Parity-mixed coupled-cluster formalism for computing parity-violating amplitudes

The field of parity violation started with the seminal paper by Lee and Yang [1] and the discovery of parity non-conservation (PNC) in nuclear $\beta$-decay [2]. Shortly after, the possibility of measuring atomic parity violation (APV) as a low-energy test for the Standard Model (SM) was investigated by Zel’dovich [3], whose consideration for hydrogen suggested that the effects were too small to be observable. The situation changed when the Bouchiats [4, 5, 6] demonstrated that APV effects scale as $Z^{3}$, where $Z$ is the nuclear charge, thus reopening the case for observing them in heavy neutral atoms. Following a proposal by Khriplovich [7], APV effects were first observed in bismuth by Barkov and Zolotorev [8]. Following this discovery, several APV experiments were performed for cesium [9, 10, 11, 12, 13], bismuth [14], lead [15, 16], thallium [17, 18] and ytterbium [19]. New APV experiments are underway or in the planning stage [20, 21, 22, 23, 24, 25, 26, 27] (see also the review [28] and references therein), with the aim of attaining a $\sim 0.1\%$ accuracy in ¹³³Cs.

The APV measurements are usually interpreted in terms of the nuclear weak charge $Q_{W}$, which is related to the measured PNC amplitude, $E_{\rm PV}$, via $E_{\rm PV}=k_{\rm PV}Q_{W}$, where $k_{\rm PV}$ is an atomic-structure factor. One wishes to compare the experimentally obtained value of $Q_{W}$ with the value predicted by the SM. For this purpose, the quantity $k_{\rm PV}$ should be known with a better accuracy than that of the amplitude $E_{\rm PV}$, thus yielding an accurate estimate of $Q_{W}$. This approach has been so far most successful in ¹³³Cs, due to its large nuclear charge, $Z=55$, and its relatively simple atomic structure with one valence electron above a closed Xe-like core [29]. Part of the success was also due to the fact that ¹³³Cs is used in the primary frequency standard, providing a wealth of information on its basic atomic properties.

In ${}^{133}\mathrm{Cs}$, the experimental uncertainty for the $6S_{1/2}-7S_{1/2}$ PNC amplitude eventually reached 0.35% [10]. The most accurate theoretical computations for the atomic-structure factor $k_{\rm PV}$ was built upon the relativistic Many-Body Perturbation Theory (MBPT), a systematic order-by-order approach which includes electron correlations. Certain classes of MBPT diagrams can be summed to all orders, taking advantage of the underlying topology of the diagrams. In the late 1980s and early 1990s, the accuracy of MBPT calculations for $k_{\rm PV}$ was estimated to be at the level of 1% [30, 31, 32, 33].

A later re-analysis [34], based on an improved theory-experiment agreement with the new measurements of atomic properties, reduced the theoretical uncertainty of Refs. [30, 31, 32, 33] to the level of 0.4%. The deduced value for the ${}^{133}\mathrm{Cs}$ weak charge differed by $2.5\sigma$ from the SM prediction, thus suggesting new physics beyond the SM [35, 36, 37, 38]. However, the inclusion of Breit [39, 40, 41] and QED radiative corrections [42, 43, 44, 45] brought the ${}^{133}\mathrm{Cs}$ result back into the essential $1\sigma$ agreement with the SM. Clearly, the answer to whether the ${}^{133}\mathrm{Cs}$ PNC result confirms the SM or hints at new physics very much depends on the quality of theoretical atomic calculation for the atomic-structure factor $k_{\rm PV}$. In the works [39, 40, 41, 42, 43, 45], the theoretical uncertainty stood at $0.5\%$, still larger than the $0.35\%$ experimental error bar.

Since the early 2000s, the theoretical error bar has been (and still is) dominated by the uncertainty of solving the basic many-body problem of atomic structure. Further progress in improving the theoretical accuracy was reported in the late 2000s [46, 47]. These calculations built upon the ab initio relativistic coupled-cluster (CC) scheme [32]. While the calculations of Refs. [31, 32, 33] were complete through the third order of MBPT for matrix elements, the scheme in Refs. [46, 47] was complete through the fourth order of MBPT. References [46, 47] have reduced the theoretical uncertainty in the ${}^{133}\mathrm{Cs}$ atomic-structure factor $k_{\rm PV}$ to 0.27%. The final value of the ${}^{133}\mathrm{Cs}$ weak charge extracted from this calculation was in an agreement with the SM prediction, placing strong constraints on a variety of new physics scenarios.

The works [32, 33, 46, 47] used a sum-over-states approach to calculate the PNC $6S_{1/2}-7S_{1/2}$ transition amplitude in ${}^{133}\mathrm{Cs}$

In this second-order expression, $\ket{nL_{J}}$ stands for various states of the ${}^{133}\mathrm{Cs}$ atom, with $n$ being the principal quantum number, $L$ the orbital angular momentum, and $J$ the total angular momentum. These are the true many-body eigen-states of the parity-proper (PP) atomic Hamiltonian. Further, $D_{z}$ is the $z$-component of the electric dipole operator ${\bf D}\equiv\sum_{i}\mathbf{d}_{i}=-\sum_{i}e{\bf r}_{i}$ and $H_{W}=\sum_{i}h_{W}(i)$ is the $P$-odd electron-nucleus weak interaction with the single-electron operator $h_{W}$ having the form

Here, $G_{F}=2.2225\times 10^{-14}$ a.u. is the Fermi constant of the weak interaction, $Q_{W}$ is the weak nuclear charge, $\rho(r)$ is the nuclear neutron density (see Ref. [48] for a discussion of neutron skin effects) and $\gamma_{5}$ is the conventional Dirac matrix.

The largest contributions to $E_{\rm PV}$ in the sum-over-states expression (1) come from terms with $n=6,7,8,9$ (the “main” contribution). In Refs. [46, 47], the required many-body states were computed using the CC approximation including singles, doubles and valence triples (CCSDvT). The computed CCSDvT wave functions were subsequently used to compute the dipole and weak interaction matrix elements entering Eq. (1). Residual contributions to Eq. (1) come from intermediate states with $n\geq 10$ (the “tail” contribution) and core-excited states. These residual contributions are sub-dominant and were evaluated using less accurate methods, having an estimated uncertainty of $10\%$.

In a later work [49], the value of the residual contributions was reevaluated and Ref. [49] claimed a contribution of core-excited states to $E_{\rm PV}$ having an opposite sign as compared to the analyses of both Refs. [32, 33] and Refs. [46, 47]. The ${}^{133}\mathrm{Cs}$ weak charge extracted from the revised atomic-structure factor is $1.5\sigma$ away from the SM value, thus relaxing Refs. [46, 47] constraints on new physics. In addition, Ref. [49] raised the theoretical uncertainty in the atomic structure factor $k_{\rm PV}$ back to 0.5%, above the experimental error bar on $E_{\rm PV}$.

The latest Dalgarno-Lewis-type coupled-cluster computations [50] support both the sign and the value of the core-excited state contributions of Refs. [46, 47]. However, as of now, a clear understanding of why the two approaches, Refs. [49] and [46, 47], lead to core-excited state contribution of opposite signs is still lacking. Furthermore, objections [51] have been raised with regards to the error estimates of Ref. [50].

It may be observed that the disagreement between Refs. [46, 47] and [49] arose due to the artificial separation into the “main” and “tail” contributions characteristic of the sum-over-states method [29, 28]. In this paper, we seek to directly include the weak interaction into the single-particle atomic Hamiltonian, thus avoiding this artificial separation, and treat all the intermediate states on equal high-precision footing. In this approach, the single-electron eigen-states of the modified Hamiltonian will already have a parity-mixed (PM) character. The MBPT calculations of the PM many-body wave functions $\ket{6S_{1/2}^{\prime}}$ and $\ket{7S_{1/2}^{\prime}}$ can be carried out in a conventional fashion using this PM singe-electron basis. This PM approach was first suggested in Ref. [52] and carried through to all second-order MBPT corrections in Ref. [53]. In this paper, we extend it to a more-complete CC method. Once the PM many-body states are computed, the PNC amplitude can be expressed simply as,

avoiding the summation over intermediate states altogether.

In addition, the lowest-order Dirac-Hartree-Fock (DHF) result in this PM approach is only 3% away from the more accurate CCSDvT value. This is to be compared with the traditional parity-proper (PP) DHF result which is off by 18%. This indicates that the correlation corrections in the PM approach are substantially smaller than in the conventional PP method. Depending on the MBPT convergence pattern, one can generically anticipate an improved theoretical accuracy.

Another important point is that in the sum-over-states approach employed in Refs. [46, 47], the theoretical uncertainty budget of $E_{\mathrm{PV}}$ included comparable contributions of the accurately computed low-lying states (in the CCSDvT approach) and of the less-accurate highly-excited and core-excited states. Our method would allow us to treat all of these contributions on the same high-accuracy CCSDvT footing, thus improving the overall theoretical uncertainty even without the potentially reduced role of the correlation corrections.

To follow this program through in the context of the CC method, one requires a numerically complete set of PM orbitals (single-particle states) $\{\psi^{\prime}_{i}\}$. Generating such PM basis sets and quantifying their numerical accuracy is one of the goals of this paper. A PM basis set has to be obtained in the modified DHF potential of the Xe-like core which includes the weak interaction (2). Considering the increased numerical accuracy demanded of the quality of basis sets, we employ the dual-kinetic-balance B-splines basis sets [54, 55] which are more numerically robust and have the correct behaviour inside the finite nucleus compared to the B-spline basis sets originally used by the Notre Dame group [56].

Once a PM basis set is obtained, one may proceed to computing matrix elements of various operators such as the one-body dipole operator $z^{\prime}_{ij}$ and the two-body inter-electron Coulomb interaction $g^{\prime}_{ijkl}$ in the new PM-DHF basis. With these computed matrix elements, the MBPT and CCSDvT expressions can be evaluated. As we will show, all the matrix elements in the PM basis can be decomposed into real and imaginary parts with opposite parities (with the conventional choice of radial wave functions being real-valued). Then all the information about opposite-parity admixtures is contained in the imaginary parts of various MBPT expressions. This greatly simplifies the formalism and only requires only minor modifications to already developed and tested MBPT codes. We demonstrate the utility of this technique for the random-phase approximation (RPA) subset of MBPT diagrams and discuss a strategy for applying these ideas in the more-complete CC calculations.

The paper is organized as follows. In Sec. II, we briefly present the set-up of our problem. Although this section does not contain new results, it serves as a starting point for our main discussion and a mean to define our notations. We will also derive, in Sec. III, the second-quantized form of the parity operator which will be useful in deriving selection rules for our PM-CC method. Section III is followed by Sec. IV, where we present several methods through which a basis of PM single-electron orbitals may be obtained. In Sec. V, we present the PM matrix elements of one- and two-body operators computed using the obtained PM single-electron orbitals. In Sec. VI, we illustrate how these PM matrix elements can be used in an RPA calculation of the PNC amplitude. The generalization to the CC method is discussed in Sec. VII. Finally, Sec. VIII draws conclusions and presents an outlook for our future work. The paper contains several appendixes which provide further technical details. Unless specified otherwise, the atomic units, $|e|=m_{e}=\hbar=1$, are used.

In this section, we lay out the theoretical framework for computing the PNC amplitude in an atom. The material presented here is not new but serves as a starting point and a mean to define our notations.

Let us begin by considering the Hamiltonian of the atomic electrons propagating in the combined PP Coulomb potential $\sum_{i}V_{\rm nuc}(r_{i})$ and the $P$-odd electron-nucleus interaction $H_{W}=\sum_{i}h_{W}(i)$. Here, $i$ labels all the atomic electrons. The full electronic Hamiltonian $H^{\prime}$ may be decomposed into

where $U^{\prime}(r_{i})$ is some single-electron potential to be specified later. We use the prime on $h^{\prime}_{0}$ and $V^{\prime}_{c}$ to distinguish them from the PP Hamiltonian $h_{0}=c\boldsymbol{\alpha}_{i}\cdot{\bf p}_{i}+m_{e}c^{2}\beta_{i}+V_{\rm nuc}(r_{i})+U(r_{i})$ and the PP $e-e$ interaction $V_{c}=\sum_{i\neq j}e^{2}/(2|{\bf r}_{i}-{\bf r}_{j}|)-\sum_{i}U(r_{i})$. For brevity, we suppressed the positive-energy projection operators for the two-electron interactions (no-pair approximation).

As usual, we assume that the energies $\varepsilon^{\prime}_{i}$ and orbitals $\psi^{\prime}_{i}$ of the unperturbed single-electron Hamiltonian $h^{\prime}_{0}$ are known. Note that since the weak interaction is a pseudo-scalar, the total angular momentum $j_{i}$ and its projection $m_{i}$ remain good quantum numbers, while the parity is no longer conserved. For example, $p_{1/2}$ and $s_{1/2}$ orbitals or $d_{5/2}$ and $f_{5/2}$ orbitals of $h_{0}$ are mixed to form eigen-states of $h^{\prime}_{0}$. The many-body eigen-states $\Psi^{\prime}$ of $H^{\prime}$ are then expanded over antisymmetrized products of the PM one-particle orbitals $\psi^{\prime}_{i}$. In MBPT, one obtains these eigen-states by treating the residual $e-e$ interaction $V^{\prime}_{c}$ as a perturbation.

As the next step, we express the terms in Eq. (4) in second quantization. Let us denote by $a^{\prime\dagger}_{i}$ and $a^{\prime}_{i}$ the creation and annihilation operators associated with the one-particle eigen-state $\psi^{\prime}_{i}$ of $h^{\prime}_{0}$. We will follow the indexing convention that core electron orbitals are denoted by the letters at the beginning of alphabet $a,b,c,\dots$, while valence electron orbitals are denoted by $v,w,\dots$, and the indices $i,j,k,\dots$ refer to an arbitrary orbital, core or excited (including valence states). The letters $m,n,p,\dots$ are reserved for those orbitals unoccupied in the core (these could be valence orbitals).

The operators $H^{\prime}_{0}\equiv\sum_{i}h^{\prime}_{0}(i)$ and $V^{\prime}_{c}$ may then be written as

where $N$ denotes normal ordering and $V^{\prime}_{\rm HF}$ is the PM-DHF potential, whose matrix elements are defined by

with $\tilde{g}^{\prime}_{ijkl}\equiv g^{\prime}_{ijkl}-g^{\prime}_{ijlk}$ being the anti-symmetrized combination of the Coulomb matrix elements,

The irrelevant constant offset energy term $\sum_{a}\left(V^{\prime}_{\rm HF}/2-U^{\prime}\right)_{aa}$ has been omitted in Eq. (5).

Notice that the choice $U^{\prime}=V^{\prime}_{\rm HF}$ causes the first term in $V^{\prime}_{c}$ in Eq. (5) to vanish, significantly reducing the number of MBPT contributions. In addition, since our final goal is to implement the CCSDvT scheme that has been originally built on DHF potential, we fix $U^{\prime}=V^{\prime}_{\rm HF}$.

With $U^{\prime}=V^{\prime}_{\rm HF}$ fixed, we now consider the correlation corrections to the independent-particle wave functions. Consider a univalent atom, e.g, ${}^{133}\mathrm{Cs}$, with a single valence electron above the closed-shell core. The zeroth-order wave function may be expressed as $\ket{\Psi^{\prime(0)}_{v}}=a^{\prime\dagger}_{v}\ket{0^{\prime}_{c}}$, where $\ket{0^{\prime}_{c}}$ represents the filled Fermi sea of the atomic core (again, the prime indicates that the single-particle orbitals are of PM character).

To the first order in the residual interaction $V^{\prime}_{c}$, the many-body correction $\delta\Psi^{\prime}_{v}$ to $\Psi^{\prime(0)}_{v}$ has the form

where we used the notation $\varepsilon^{\prime}_{ij}\equiv\varepsilon^{\prime}_{i}+\varepsilon^{\prime}_{j}$. Here, the first term describes a valence electron being promoted to an excited state orbital $m$ with a simultaneous particle-hole excitation of the core (this is so-called valence double excitation, $D_{v}$, in the language of the CC method). The second contribution is a double particle-hole excitation from the core with valence electron being a spectator (core double excitation, $D_{c}$). The anti-commutation relations for creation and annihilation operators assure that the electrons in the second term do not get excited into the valence orbital, i.e. the Pauli exclusion principle is built into the formalism automatically.

With $\ket{\delta\Psi^{\prime}_{v}}$, one can compute the second-order correction to the matrix element of a one-electron operator $T=\sum_{ij}t^{\prime}_{ij}a^{\prime\dagger}_{i}a^{\prime}_{j}$. Once again, the primed quantities refer to the PM orbitals used in computing the matrix elements $t^{\prime}_{ij}\equiv\langle i^{\prime}|t|j^{\prime}\rangle$. The operator $T$ can be, for example, the electric dipole operator $\bf D$. The correction to the matrix element between two valence many-body states $\ket{\Psi^{\prime}_{w}}$ and $\ket{\Psi^{\prime}_{v}}$, $w\neq v$, has the form

where $\omega\equiv\varepsilon^{\prime}_{w}-\varepsilon^{\prime}_{v}$ and $\tilde{g}^{\prime}_{ijkl}\equiv g^{\prime}_{ijkl}-g^{\prime}_{ijlk}$. Equations (II) and (II) are, of course, well-known results [57] with the only difference of using the mixed-parity basis instead of the conventional PP basis.

In Eq. (II), we sum over the core orbitals $a$ and the excited orbitals $n$. Each orbital $\psi^{\prime}_{i}$ is characterized by a principal quantum number $n_{i}$, a total angular momentum $j_{i}$, and its projection $m_{i}$. The sums over the magnetic quantum numbers $m_{i}$ can be carried out analytically using the rules of Racah algebra. Although the sums over $j_{i}$ are infinite, they are restricted by angular momentum selection rules which radically reduce the number of surviving terms. Moreover, the sums over total angular momenta converge well and in practice, it suffices to sum over a few lowest values of $j_{i}$. The sums over the principal quantum numbers $n_{i}$ involve, on the other hand, summing over the infinite discrete spectrum and integrating over the continuum. In the basis set method, these infinite summations are replaced by summations over a finite-size pseudo-spectrum [58, 59, 54, 55].

The basis orbitals in the pseudo-spectrum are obtained by placing the atom in a sufficiently large cavity and imposing boundary conditions at the cavity wall and at the origin (see Ref. [55] for further details on a dual-kinetic-basis B-spline sets used in our paper). For each value of $j_{i}$, one then finds a discrete set of $2N$ orbitals, $N$ from the Dirac sea and the remaining $N$ with energies above the Dirac sea threshold (conventionally referred to as “negative” and “positive” energy parts of the spectrum in analogy with free-fermion solutions).

If the size of the cavity is large enough, typically about $40a_{0}/Z$ where $a_{0}$ is the Bohr radius, the low-lying orbitals with positive energies map with a good accuracy to the discrete orbitals of the exact DHF spectrum. Higher-energy orbitals do not closely match their physical counterparts. Nevertheless, since the pseudo-spectrum is complete, it forms a basis set for the function space spanning the cavity and thus can be used instead of the real spectrum to evaluate correlation corrections to states confined to the cavity. From now on, all single-particle orbitals $\psi^{\prime}_{i}$ are understood to be members of the B-spline basis set.

To reiterate, the parity-mixed (PM) formalism presented so far is essentially the same as in the conventional parity-proper (PP) MBPT. The only difference is that all the quantities are defined with respect to the PM orbitals $\psi^{\prime}_{i}$ instead of the PP ones. Since PM orbitals are eigen-states of total angular momentum, one can directly use the results of angular reduction for various MBPT expressions, and the existing MBPT codes require minor changes, mostly related to modifying parity selection rules and the use of Coulomb integrals in the PM basis. In Sec. III, we derive the second-quantized form of the parity operator which will be useful for deriving the PM selection rules and in Sec. IV, we present several methods through which the PM orbitals may be generated in practice.

Since the MBPT derivations are built on the second-quantization formalism, in this section we derive the second-quantized form of the parity operator $\Pi$ to be used in deriving parity selection rules (see Appendix B). Parity transformation is defined by $\mathbf{r}_{i}\rightarrow-\mathbf{r}_{i}$ for all the $N_{e}$ electrons in the system. Consider a PP state (Slater determinant) $\ket{\Psi^{m_{1}\ldots m_{\mu}}_{a_{1}\ldots a_{\mu}}}$ composed of orbitals of definite parity. This many-body state is obtained by removing $\mu=0,\ldots,N_{e}$ electrons $a_{1},\ldots,a_{\mu}$ from the reference state $\ket{\Psi^{(0)}_{v}}$ while adding the same number of excited electrons $m_{1},\ldots,m_{\mu}$. Notice that in this notation the valence orbital is treated as initially occupied and, thereby, $v$ can be one of the labels $a_{1},\ldots,a_{\mu}$. In the second quantization,

We emphasize that the creation and annihilation operators in Eq. (10) are the PP ones.

Since the PP Hamiltonian $H_{0}\equiv\sum_{i}h_{0}(i)$ is invariant under spatial reflection, it commutes with the parity operator,$[H_{0},\Pi]=0$. As a result, the states $\ket{\Psi_{v}^{(0)}}$ and $\ket{\Psi^{m_{1}\ldots m_{\mu}}_{a_{1}\ldots a_{\mu}}}$, being eigen-states of $H_{0}$, are also an eigen-states of the parity operator $\Pi$. Furthermore, since $\ket{\Psi_{v}^{(0)}}$ and $\ket{\Psi^{m_{1}\ldots m_{\mu}}_{a_{1}\ldots a_{\mu}}}$ are antisymmetrized products of single-electron orbitals, their eigen-values with respect to $\Pi$ equal the products of the parities of their constituents

where we have used the fact that the closed-shell core has even parity.

To transform an operator into the second-quantized form, we recall the conventional formula

Its proof relies on the identity resolution (closure relation) for a complete orthonormal basis $I=\sum_{\alpha}\ket{\alpha}\bra{\alpha}=\sum_{\beta}\ket{\beta}\bra{\beta}$. For a system of identical particles, however, one needs to proceed with caution due to the possibility of permutations of orbitals in $\ket{\Psi^{m_{1}\ldots m_{\mu}}_{a_{1}\ldots a_{\mu}}}$. Indeed, in the many-fermion case, the orthonormality condition reads

where the generalized Kronecker delta is defined as

For many-fermion systems, the general closure relation is given in Ref. [60]. Since $\Pi$ is a diagonal operator, the general identity resolution [60] simplifies to

where $\{a\}$ and $\{m\}$ denote strings of orbital labels in $\ket{\Psi^{m_{1}\ldots m_{\mu}}_{a_{1}\ldots a_{\mu}}}$.

Sandwiching $\Pi$ in between two identity operators and using the closure relation (15), we find

Using the eigenvalue equation (11b), we obtain the following result for the matrix element of $\Pi$

where we used the orthonormality relation (13). Substituting Eq. (14) into Eq. (III), we finally obtain

In this section, we demonstrate how the PM single-particle wave functions $\psi^{\prime}_{i}\equiv\psi^{\prime}_{n_{i}j_{i}m_{i}}$ may be obtained. They are the solutions to the PM-DHF equation

Here, $n_{i}$ is the principal quantum number, $j_{i}$ is the total angular momentum and $m_{i}$ is the projection of $j_{i}$ on a quantization axis.

Our goal is to expand $\psi^{\prime}_{i}$ in terms of the PP orbitals $\psi^{P}_{i}\equiv\psi^{P}_{n_{i}\ell_{i}j_{i}m_{i}}$, which are solutions to the conventional DHF equation

Note that besides the principal quantum number $n_{i}$, the total angular momentum $j_{i}$, and the magnetic quantum number $m_{i}$, we have characterized the PP orbital $\psi^{P}_{i}$ with an extra quantum number, the orbital angular momentum $\ell_{i}$, which indicates that $\psi^{P}_{i}$ has a definite parity $(-1)^{\ell_{i}}$.

The two DHF potentials $V^{\prime}_{\mathrm{HF}}$ and $V_{\mathrm{HF}}$ depend on the core orbitals. Since core orbitals are self-consistent solutions of these DHF equations, the PM and PP core orbitals differ. Therefore, the effects of the weak interaction on the single-electron orbitals are contained in the difference of the two Hamiltonians,

which is a pseudo-scalar interaction, preserving rotational symmetry but spoiling mirror symmetry (this is why we have used the quantum numbers $n_{i}$, $j_{i}$ and $m_{i}$ but not $\ell_{i}$ to index the PM orbital $\psi^{\prime}_{i}$). This suggests the following parameterization [61] for the solutions to the PM-DHF equations (19),

where

is a dimensionless factor characteristic of the strength of the weak interaction and its numerical value is given for ${}^{133}\mathrm{Cs}$. In Eqs. (22), $\kappa_{i}=\left(j_{i}+\frac{1}{2}\right)(-1)^{j_{i}+\ell_{i}+1/2}$ is a relativistic angular quantum number that encodes the values of both the total and orbital angular momenta $j_{i}$ and $\ell_{i}$. From the definition of the relativistic angular quantum number, flipping the parity $(-1)^{\ell_{i}}$ of the orbital while preserving the total angular momentum $j_{i}$ is equivalent to changing $\kappa_{i}\rightarrow-\kappa_{i}$, as presented in the parameterization of $\bar{\psi}_{i}$, Eq. (22g).

It is appropriate to pause here and introduce a point of semantics. Although the orbital $\psi^{\prime}_{i}$ does not have a definite parity, one can nevertheless speak of its “nominal parity”, defined as that of the component $\psi_{i}$, which is not suppressed by the factor $\eta$. In the light of Eq. (22a), we shall refer to $\psi_{i}$ as the “real” component and $\bar{\psi}_{i}$ as the “imaginary” component of $\psi^{\prime}_{i}$. In what follows, in particular when discussing MBPT and the CC formalism, we shall refer to the nominal parity of a PM single-electron orbital, meaning that of its real component. The nominal parity of $\psi^{\prime}_{i}$ is thus $(-1)^{\ell_{i}}$.

The combination (22a) clearly demonstrates the admixing of the opposite-parity orbital $\bar{\psi}_{i}$ with $\kappa_{i}\rightarrow-\kappa_{i}$ to the reference orbital $\psi_{i}$. With the imaginary unity factored out in the admixture component $\bar{\psi}_{i}$ and with the conventional definition of the spherical spinors $\Omega_{\kappa m}$ [62], all the radial wave functions $P_{n_{i}\kappa_{i}}$, $Q_{n_{i}\kappa_{i}}$, $\bar{P}_{n_{i}\kappa_{i}}$, and $\bar{Q}_{n_{i}\kappa_{i}}$ can be chosen to be real-valued. The rest of this section will be devoted to finding these radial components. In what follows, we will assume a parameteriziton for the PP solutions to Eq. (LABEL:Normal_DHF) similar to that presented in Eqs. (22d), namely

Where there is no risk of confusion, we will abbreviate $P_{n_{i}\kappa_{i}}$ ($Q_{n_{i}\kappa_{i}}$) to $P_{i}$ ($Q_{i}$) and $P^{P}_{n_{i}\kappa_{i}}$ ($Q^{P}_{n_{i}\kappa_{i}}$) to $P^{P}_{i}$ ($Q^{P}_{i}$). The energy eigenvalues $\epsilon^{\prime}_{n_{i}j_{i}}$ and $\epsilon_{n_{i}l_{i}j_{i}}$ in Eqs. (19) and (LABEL:Normal_DHF) will also be abbreviated to $\epsilon^{\prime}_{i}$ and $\epsilon_{i}$, respectively.