Comment on “New physics constraints from atomic parity violation in ¹³³Cs”

QED radiative corrections in the strong Coulomb field of the nucleus make a significant contribution to $E_{\rm PV}$, $\lesssim$ 1%. These have been calculated before Johnson et al. (2001); Dzuba et al. (2002); Milstein et al. (2002); *Milstein2002a; *Milstein2003; Kuchiev and Flambaum (2002); Kuchiev (2002); *KuchievQED03; Sapirstein et al. (2003); Shabaev et al. (2005a, b); Flambaum and Ginges (2005); Roberts et al. (2013a) and are well established. It is said in the Letter Sahoo et al. (2021) that one of the key improvements is the treatment of these QED corrections. However, details of the QED calculation are not presented in the Letter, and the reader is directed to the unpublished manuscript Sahoo and Das (2020) for explanation ¹¹1Note that the reference to Sahoo and Das (2020) in the Letter Sahoo et al. (2021) is incorrect, linking to an unrelated arXiv paper. There it is said that the self-energy QED correction to $E_{\rm PV}$ (and to other atomic properties) is accounted for by including the radiative potential Flambaum and Ginges (2005); Ginges and Berengut (2016a) into the Hamiltonian from the start ²²2Vacuum polarization is included using the standard Uehling potential, and a simplified form of the Wichmann-Kroll potential from Ref. Dzuba et al. (2002); Flambaum and Ginges (2005)., which the authors claim to be a more rigorous approach compared to previous calculations.

The radiative potential method Flambaum and Ginges (2005) enables the accurate inclusion of self-energy corrections to the energies and wavefunctions of many-electron atoms. It may also be used to account for QED corrections to matrix elements of external fields whose operators act at radial distances much larger than the electron Compton wavelength, $r\gg e\hbar/(m_{e}c)$, e.g., the E1 field. However, this is not the case for operators that act at small distances, including the weak and hyperfine interactions. We illustrate this in Fig. 1. For the E1 interaction, the dominant contribution is given by the left and right diagrams, which may be accounted for by using the radiative potential method. However, for the weak and hyperfine interactions, other contributions are important. In particular, the middle vertex diagram – where the external field is locked inside the photon loop – simply cannot be accounted for using this method. We refer the reader to the original Flambaum and Ginges (2005) and subsequent Roberts et al. (2013a) works for details on the applicability of the radiative potential method.

The QED correction to the full Cs APV amplitude (involving both E1 and weak interactions) was determined in Refs. Flambaum and Ginges (2005); Shabaev et al. (2005a). In Ref. Flambaum and Ginges (2005), the radiative potential method was used to calculate corrections to the E1 matrix elements and energy denominators in the sum (1), with QED corrections to weak matrix elements $\langle s|h_{w}|p_{1/2}\rangle$ taken from previous works Milstein et al. (2002); *Milstein2002a; *Milstein2003; Kuchiev and Flambaum (2002); Kuchiev (2002); *KuchievQED03; Sapirstein et al. (2003). In Ref. Shabaev et al. (2005a), Shabaev et al. calculated the total correction by applying a rigorous QED formalism. The results of Refs. Shabaev et al. (2005a) and Flambaum and Ginges (2005) are in excellent agreement, $-0.27(3)$% and $-0.32(3)$%, respectively ³³3The QED results of both Refs. Shabaev et al. (2005a) and Flambaum and Ginges (2005) were misattributed in Table III of the Letter Sahoo et al. (2021); these papers were not cited in Ref. Sahoo et al. (2021)..

It is unclear how the authors of Sahoo et al. (2021) arrive at a QED correction of $-0.4\%$ for the weak matrix elements and $-0.3\%$ for $E_{\rm PV}$, in agreement with existing calculations Milstein et al. (2002); *Milstein2002a; *Milstein2003; Kuchiev and Flambaum (2002); Kuchiev (2002); *KuchievQED03; Sapirstein et al. (2003); Flambaum and Ginges (2005); Shabaev et al. (2005a, b); Roberts et al. (2013a), given the important short-range effects, including the vertex contribution, have been omitted. In an attempt to reproduce the results of Ref. Sahoo et al. (2021), we calculate the radiative potential value for the QED correction to weak matrix elements, including vacuum polarization. The result is $-2.1\%$, too large by a factor of five compared to the correct calculations, demonstrating the importance of the missed short-range effects. This difference amounts to a change in $E_{\rm PV}$ that is nearly six times the atomic theory uncertainty claimed in Ref. Sahoo et al. (2021).

In the Letter Sahoo et al. (2021), calculations of hyperfine constants are performed to test the accuracy of the wavefunctions in the nuclear region, crucial for assessing the accuracy of APV calculations (see Refs. Ginges et al. (2017); Ginges and Volotka (2018); Roberts and Ginges (2020, 2021) for recent studies of the nuclear magnetization distribution for Cs). By demonstrating excellent agreement with experiment, the authors conclude the accuracy of their wavefunctions is high, and so estimate a tremendously small uncertainty for the APV calculation. However, it appears that serious omissions have been made in the hyperfine calculations.

As for $E_{\rm PV}$, the vertex and short-range contributions to QED corrections to hyperfine constants are important Sapirstein and Cheng (2003); Ginges et al. (2017) (see also Blundell et al. (1997); Sunnergren et al. (1998); Artemyev et al. (2001); Volotka et al. (2008, 2012)). Moreover, the magnetic loop vacuum polarization correction also gives a significant contribution Sapirstein and Cheng (2003); Ginges et al. (2017). In the Letter Sahoo et al. (2021), the radiative potential method is employed, with no account for these contributions. Given this, it is unclear how the authors of Sahoo et al. (2021); Sahoo and Das (2020) arrive at a correction of $-0.3\%$ to the hyperfine constants for $s$ states of Cs, in good agreement with existing calculations Sapirstein and Cheng (2003); Ginges et al. (2017). To investigate this result, we again use the radiative potential method and find it gives a correction of $-1.2\%$, three times too large compared to rigorous QED calculations Sapirstein and Cheng (2003); Ginges et al. (2017), confirming the importance of the omitted effects. This difference amounts to two times the uncertainty of the hyperfine calculations (0.4%) claimed in the Letter Sahoo et al. (2021).

The contribution to $E_{\rm PV}$ coming from the (occupied) $n$ = 2–5 terms in Eq. (1) is called the “core” (or autoionization) contribution. In the Letter Sahoo et al. (2021), it is said that the main difference in the $E_{\rm PV}$ result compared to the previous calculation of Dzuba et al. Dzuba et al. (2012) stems from the opposite sign of the core contribution. The difference in core contribution between Refs. Sahoo et al. (2021) and Dzuba et al. (2012) is larger than the theoretical uncertainty claimed in the Letter Sahoo et al. (2021) and should be investigated thoroughly.

In Ref. Dzuba et al. (2012), Dzuba et al. showed that many-body effects (core polarization and correlations) have a significant impact on the core contribution, changing its sign compared to the lowest-order Hartree-Fock value; see also Ref. Roberts et al. (2015). The authors of Ref. Sahoo et al. (2021) claim their result confirms the core calculation of Ref. Porsev et al. (2009); *Porsev2010 and agrees with the result of Ref. Blundell et al. (1990); *Blundell1992. However, in both of those works, the core contribution was evaluated in the lowest-order approximation.

Here, we re-examine the core contribution in detail in an attempt to elucidate the source of this discrepancy. We include core polarization using the time-dependent Hartree-Fock (TDHF) method Dzuba et al. (1984), in which the single-particle operators are modified: $d_{z}\to{\tilde{d}_{z}}=d_{z}+\delta V_{d}$, and ${h_{w}}\to{\tilde{h}_{w}}={h}_{w}+\delta V_{w}$. The $\delta V$ corrections are found by solving the set of TDHF equations for all electrons in the core Dzuba et al. (1984). We obtain the corrections to lowest-order in the Coulomb interaction by solving the set of equations once, and to all-orders by iterating the equations until self-consistency is reached Dzuba et al. (1984) (equivalent to the random-phase approximation with exchange, RPA Johnson et al. (1980)). The equations for $\delta V_{d}$ are solved at the frequency of the $6S$–$7S$ transition (see Roberts et al. (2013b) for a numerical study). We account for correlation corrections using the second-order Dzuba et al. (1987) and all-orders Dzuba et al. (1988); *DzubaCPM1989plaEn; *DzubaCPM1989plaE1 correlation potential methods (see also Dzuba et al. (2002)).

The core contribution arises as the sum of two terms, due to the weak-perturbation of $6s$ and $7s$ states, respectively. These have similar magnitude though opposite sign, and strongly cancel, meaning numerical error may be significant. We test the numerical accuracy in a number of ways. Firstly, we vary the number of radial grid points used for solving the differential equations, and vary the number of basis states used in any expansions. We find numerical errors stemming from grid/basis choices can easily be made insignificant. More importantly, we have three physically equivalent, but numerically distinct, ways to compute $E_{\rm PV}$:

where $\delta\psi$ and $\Delta\psi$ are corrections to the valence wavefunctions ($\psi$) due to the time-independent weak interaction, and the time-dependent E1 interaction, respectively. These are called the sum-over-states (2), weak-mixed-states (3), and E1-mixed-states (4) methods ⁴⁴4These formulas exclude the double-core-polarization effect, which is very small for Cs, and has been studied in detail in Ref. Roberts et al. (2013b). The sign change of the core cannot be explained by the double-core-polarization correction, which even if entirely assigned to the core contribution is a factor of two too small to account for the difference Roberts et al. (2013b)..

In the sum-over-states method, a B-spline basis (e.g., Johnson et al. (1988); Beloy and Derevianko (2008)) is used to sum over the set of intermediate states. In contrast, the mixed-states approach does not require a basis at all; the $\delta$ and $\Delta$ corrections are found by solving the differential equations Dalgarno and Lewis (1955):

where $h$ is the single-particle atomic Hamiltonian, and $\nu$ is the $6S$–$7S$ transition frequency. In the mixed-states approach, the core contribution is found by projecting the corrections $\delta\psi$ and $\Delta\psi$ onto the core states, while in the sum-over-states method it is found by restricting the sum to include only core states. Note that the numerics involved in solving each of the above equations is significantly different, and the coincidence of results is indicative of high numerical accuracy. Even with moderate choice for the radial grid, we find the results of the two mixed-states methods agree to parts in $10^{8}$, and the mixed-states and sum-over-states methods agree to parts in $10^{7}$, demonstrating excellent numerical precision and completeness of the basis.

Our calculations of the core term are summarized in Table 1. The sign change in the core contribution is mostly due to polarization of the core by the external E1 field. This is sensitive to the frequency of the E1 field. While correlations beyond core polarization are important, they affect both terms in roughly the same manner; the core term and its sign are robust to the treatment of correlations. We also performed calculations for the $7S$-$6D_{3/2}$ $E_{\rm PV}$ for ²²³Ra⁺ to test against previous calculations; at the RPA level, we find the core contribution to be 6.81 [in units $-10^{-11}i(-Q_{w}/N)\,|e|a_{B}$], in excellent agreement with the result 6.83 of Ref. Pal et al. (2009) (see also Dzuba et al. (2001b); Wansbeek et al. (2008)). It is unclear why the sign of the result of Ref. Sahoo et al. (2021) remains the same as the Hartree-Fock value, however, we note that it may not be straight forward to compare individual contributions across different methods as discussed in Refs. Wieman and Derevianko (2019); Safronova et al. (2018).

For the above reasons, we are not convinced the result presented in the Letter Sahoo et al. (2021) is an improved value for the Cs $E_{\rm PV}$. We conclude that the most reliable and accurate values that have been obtained to date are: $E_{\rm PV}=0.898(5)$ Dzuba et al. (2002); Flambaum and Ginges (2005) and $E_{\rm PV}=0.8977(40)$ Porsev et al. (2009); Dzuba et al. (2012), in units $-10^{-11}i(-Q_{w}/N)\,|e|a_{B}$, which agree precisely and were obtained using different approaches. These results are also in excellent agreement with previous calculations Dzuba et al. (1989c); Blundell et al. (1990); *Blundell1992; Kozlov et al. (2001); Shabaev et al. (2005b), though in disagreement with the result of the Letter Sahoo et al. (2021).

Acknowledgments— We thank V. A. Dzuba and V. V. Flambaum for useful discussions. This work was supported by the Australian Government through ARC DECRA Fellowship DE210101026 and ARC Future Fellowship FT170100452.