Relativistic effects in search for new intra-atomic force with isotope shifts

Isotope shift of atomic spectra is considered as a probe of new interaction between electrons and neutrons in atoms. We employ the method of seeking a breakdown of King's linearity in the isotope shifts of two atomic transitions. In the present work, we evaluate the magnitudes of the nonlinearity using relativistic wave functions and the result is compared with that of nonrelativistic wave functions in our previous work. It turns out that the nonrelativistic calculation underestimates the nonlinearity owing to the new interaction in the mass range of the mediator greater than 1 MeV. Further, we find that the nonlinearity within the standard model of particle physics is significantly magnified by the relativistic effect in the $\text{p}_{1/2}$ state. To get rid of this obstacle in the new physics search, we suggest to avoid $\text{p}_{1/2}$, and use $\text{p}_{3/2}$ instead for example.


Introduction
In developing the next-generation frequency standard, the precision of atomic spectroscopy has been considerably improved. For instance, the relative uncertainty of O(10 −18 ) is realized in a clock transition of Yb + [1]. The corresponding accuracy of frequency (energy) is of the order of mHz (10 −18 eV).
We consider a possibility to constrain or probe a new physics beyond the standard model of particle physics with such high precision spectroscopy. It is practically impossible to calculate transition frequencies of multielectron atomic systems at the level of the above accuracy. To overcome this difficulty, a method to use isotope shifts (IS) was proposed [2,3].
The typical magnitude of observed IS is the order of GHz (10 −6 eV), much smaller than optical transition energies, and such quantities are evaluated by perturbation. If nuclei were infinitely heavy and point-like, the electrons in an atom could not sense the difference of the nuclei among isotopes and no IS would exist. Thus, it is customary to decompose observed isotope shifts into the mass shift (MS), which is due to the finite nuclear mass, and the field shift (FS), due to the finite nuclear size [4]. In other words, the nuclear recoil causes the MS and the change of the nuclear Coulomb field does the FS.
The isotope shift of a transition (labeled by t) between an isotope pair A and A , ν t,A A := ν t,A − ν t,A , in the leading order perturbation is expressed by where µ A A := µ A − µ A represents the difference of the reduced masses, r 2 A A := r 2 A − r 2 A denotes that of the mean-squared radii of nuclei, and K t and F t are electronic factors. The electronic factors solely depend on the transition and are common for isotopes.
The key idea to explore new physics using IS is King's linearity [5] and its breakdown by a new interaction [2]. Suppose isotope shifts of several isotope pairs in two distinct transitions (t = 1, 2) of an atom (or an ion). It is convenient to introduce the modified isotope shift (mIS)ν t,A A := ν t,A A /µ A A so that Eq. (1) is written as We note that atomic weights are measured with typical precision of 10 −8 or better [6]. Eliminating the nuclear factor r 2 A A /µ A A , we obtain a linear relation between the mIS's of two transitions,ν where This linear relation is broken by a new force of particle exchange between an electron and a neutron. In the presence of such particle exchange, a particle shift (PS), a new source of IS, is added to the IS expression in Eq. (1), The linear relation in Eq. (3) is also modified as where and M denotes the atomic mass unit. We use ε PS to quantify the nonlinearity due to PS for a pair of transitions. We note that values of ε PS for different pairs of transitions can not be directly compared in general. In the standard model, as pointed out in Ref. [2], Z or Higgs boson leads to such a nonlinearity. However, the nonlinearities by such heavy particles are strongly suppressed and its measurement is beyond the reach of the present technology [2,7]. If there exists a light boson that couples to the electron and the neutron, the nonlinearity may be detectable depending on its mass and couplings.
One of the other conceivable sources of nonlinearity within the standard model is the FS other than the leading contribution shown in Eq. (1), e.g. a term proportional to r 4 A A [8] as mentioned in Sec. 4. Such subleading contributions are effectively included in Eq. (5) by replacing ε PS with ε PS + ε FS , where ε FS represents the FS nonlinearity. The FS nonlinearity is a potential background in the search of new physics with the IS nonlinearity, though it can be removed by using three or more transitions in principle [7].
In this work, we evaluate the PS and FS nonlinearities employing relativistic wave functions and compare the result with that of our previous work with nonrelativistic (NR) wave functions [7]. We find that the sensitivity to a new electron-neutron interaction is enhanced in the mediator mass range of about 1 to 10 MeV in our relativistic evaluation but the FS nonlinearity significantly increases as far as a p 1/2 state is involved in the relevant transitions.
The rest of the paper is organized as follows. We present the formulation to evaluate PS and FS in Sec. 2. In Sec. 3, our method to obtain relativistic wave functions is described. We give formulas to estimate PS and FS nonlinearities with the obtained wave functions in Sec. 4. The current experimental status and future prospects are given in Sec. 5. Section 6 is devoted to our conclusion. We use the natural units, = c = 1, unless otherwise stated.

Formulation of isotope shift
As in the previous work, we employ the single-electron approximation to evaluate PS and FS. One valence electron is supposed to change its state (orbital) in a transition while the state of the other electrons is kept intact. We denote the radial densities of the initial and final states of this electron by σ i (r) and σ f (r) respectively. We account for them in Sec. 3.
A new force between an electron and a neutron by exchanging a scalar or vector boson of mass m is described by the potential, where s is the mediator spin, and g n and g e denote the neutron and electron coupling constants respectively. The PS electronic factor in Eq. (4) is given by where σ t (r) := σ i (r)−σ f (r) and the radial densities are normalized as dr r 2 σ i,f (r) = 1. We note that the point-like nucleus is assumed in Eq. (7) so that it is valid for the inverse of the mediator mass much larger than the nuclear radius ∼ (100 MeV) −1 .
We also evaluate the FS with the same electron densities to estimate the magnitude of the FS nonlinearity. The FS is the variation of the transition energy owing to the difference of the Coulomb potential of two nuclei A and A . Assuming spherically symmetric nuclear charge distributions, the potential difference between isotopes is also spherically symmetric and represented by where Z is the nuclear charge, ρ A A (r) := ρ A (r)−ρ A (r), and the spherical nuclear charge distribution of an isotope A denoted by ρ A (r) is normalized as 4π dr r 2 ρ A (r) = 1. Then, the FS of transition t is given by

Relativistic wave functions
In order to evaluate the radial electron density σ i,f (r) in a relativistic manner, we solve the Dirac equation with an effective potential that represents the potential by the nucleus and the electrons other than the one involved in a transition.

Nuclear charge distribution
Since the nuclear potential difference V A A (r) vanishes outside the nuclei, the FS in Eq. (10) is governed by the wave functions inside the nuclei. Hence it is important to evaluate the wave functions near the origin taking account of the finite nuclear size. We employ the Helm distribution [9], which is the gaussian-smeared uniform sphere, The nuclear Coulomb potential of the above Helm charge distribution is given by 1 where the error function is defined by so that Erf(+∞) = 1. A more detailed description of the Helm distribution is found in Appendix B in Ref. [7]. Following Ref.

Effective potential
The electronic states in the single-electron approximation are determined by an effective potential that describes the nuclear field and the field of other electrons kept unchanged in a transition. As in our previous work, we make use of the Thomas-Fermi (TF) model, in which the ensemble of atomic electrons contributing to the effective potential is treated as a free Fermi gas in a slowly varying external potential. It gives the following effective potential for a point-like nucleus, where n is the total charge (in the unit of |e|) of the system described by the effective potential, r 0 = bx 0 with x 0 the zero of the TF function χ(x), b = [(3π) 2 /(2 7 Z)] 1/3 a B and a B := 1/αm e is the Bohr radius. The TF function satisfies with χ(0) = 1. The second boundary condition that uniquely specifies the solution is x 0 χ (x 0 ) = −n/Z with x 0 being the solution of χ(x 0 ) = 0 for a positive ion (n ≥ 1).
In the following part of the paper, we examine IS's of singly charged positive ions so that the single-electron orbitals are determined by the effective potential of n = 2. See e.g. Ref. [11] for details of the TF model. The potential in Eq. (15) is obtained for the case of the point-like nucleus. As stressed above, the finite nuclear size should not be overlooked in the FS calculation. Accordingly we modify the TF potential by subtracting the Coulomb potential of the point-like nucleus, V c (r) = −Zα/r, and adding the potential of the Helm distribution V A (r), This modification significantly alters the behavior of some wave functions inside the nucleus though it affects the spectra little.

Dirac equation
The relativistic wave function of an electron bounded in the modified TF potential in Eq. (17) is obtained by solving the eigenvalue problem of the Dirac equation, where m e is the electron mass, and α and β denote the Dirac matrices. It is well known that the eigenfunction of this central-force problem is expressed in the form of separation of variables in the spherical coordinates as where Y j 3 j (θ, ϕ) represents the spinor spherical harmonics, and G(r) and F (r) are radial wave functions. See e.g. Ref. [12].
We solve the radial equations in Eqs. (20) and (21) with the numerical method of shooting described in Ref. [13]. state is plotted in the atomic units (solid red). The NR 6p density (black dotted) is also presented for comparison. We show 6p 3/2 (red dashed) as well for a later discussion. The inset is a close look near the nucleus, whose size is saliently different from the NR 6p state, As discussed below, this difference results in an enhancement of the FS nonlinearity in the relativistic calculation. We note that the 6p 3/2 behaves in a similar way as the NR state and its possible benefit is discussed in Sec. 5.

Nonlinearity formulas
The leading contribution of the FS, F t in Eq. (1), is identified by expanding σ t (r) as σ t (r) = σ t (0) + σ t (0)r 2 /2 + · · ·. We note that, as illustrated in the inset of Fig. 1, no linear term appears provided that the nuclear charge distribution has no cusp at the origin [8].
As for the PS, we find and that the leading σ t (0) term in the case of heavy mediator disappears in the expression of the PS nonlinearity in Eq. (6) so that ε PS ∼ O(1/m 4 ) as m → ∞ [7]. In our analysis below, we numerically evaluate the integration in X t without expanding σ t (r). The FS is also obtained by numerical integration in Eq. (10) (or (22)) without expansion. We decompose the FS into the leading and subleading contributions as ). Then, the FS nonlinearity is given by as in the same manner as the PS nonlinearity. Strictly speaking, G t evaluated by subtracting the leading FS from the total FS depends on the isotope pair (A, A ). However, our numerical result shows that the dependence is weak and less than a few per cent among the isotope pairs analyzed below.

Current status and future prospects
We analyze the same experimental data of the singly charged calcium and ytterbium ions as our previous work [7], in which NR wave functions are employed. The IS in the transitions of 2 S 1/2 -2 P 1/2 (397 nm) and 2 D 3/2 -2 P  [7].
It turns out that the data satisfy the linearity within the error. Comparing the data with the IS formula including the nonlinearity, Eq. (5) with the PS nonlinearity parameter ε PS replaced by the experimental nonlinearity parameter ε exp , we obtain quantitative bounds on the nonlinearity as presented in Table 1. We also present the expected sensitivity of future experiments in Table 1 assuming the experimental precision of 1 Hz. We note that IS measurements of 1 Hz precision are conceivable since the IS of the Ca + 2 S 1/2 -2 D 5/2 (729 nm) transition is measured with a precision better than 10 Hz [18].
Our prediction of the PS nonlinearity ε PS indicated in Table 1 is that for a representative set of mediator mass and couplings in Eq. (7), m = 1 keV and |g n g e | = 1 × 10 −13 . The PS nonlinearity of Ca + for this set of parameters is the same order as the sensitivity expected in experiments of 1 Hz precision. For Yb + , the sensitivity of 1 Hz experiments is sufficiently high to probe the PS nonlinearity for the same set of parameters, provided that the FS nonlinearity is suppressed. Figure 2 shows the current constraints on the mass and couplings as well as those expected with experimental data of 1 Hz precision. The region above the upper lines are excluded by the present data, while the lower lines express the expected sensitivity of future experiments. The red (blue) lines are the case of Ca + (Yb + ). The results of the relativistic calculation in the present work are represented by the solid lines and those of the NR calculation in our previous work [7] are shown by the dash-dotted lines for comparison. We observe that the PS nonlinearity is enhanced and the sensitivity is improved significantly for m 1 MeV in the relativistic calculation of Yb + compared to the NR case. This is due to the inward attraction of relativistic wave functions. The relativistic effect is less notable for Ca + in the depicted mass range.
The peak structure that implies the loss of sensitivity at the corresponding mediator mass is due to the cancellation in the PS electronic factor in Eq. (8). Since the peak position depends on the involved wave functions, it is useful for evading the sensitivity loss to combine two or more linearity tests, e.g. Ca + and Yb + as shown in Fig. 2.
The current and expected constraints of IS nonlinearity are compared with those by other experiments and observations for each case of a scalar and vector mediator in Fig. 3. The constraints by the IS are identical for both cases. The grey shaded regions of smaller mass in both panels are excluded by the fifth force search [19,20]. The shaded regions bounded from below by the orange lines indicate the combined excluded regions by the electron g − 2 measurement [21], the dark photon search at Babar [22], the beam dump experiments [23,24] and the neutron scattering experiments [25,26,27]. The combined terrestrial bounds are not identical for the scalar and vector mediators. The left side regions of the brown lines are constrained by the stellar cooling [28]. We note that there is an uncertainty that may invalidate these stellar constraints in the relatively strong coupling regions above the brown dotted lines [29]. The black vertical bar in the right panel indicates the coupling range suggested by the 17 MeV Atomki anomaly in the 8 Be internal conversion [30,31,32]. We observe that the future IS experiments, especially Yb + , have potential sensitivity to probe the parameter space of mass range from O(0.1) to O(100) keV that has not been excluded by other experiments nor observations.
In the fifth column of Table 1, we present the FS nonlinearity parameter in Eq. (25) evaluated with the relativistic wave functions. It turns out that the FS nonlinearity  is enhanced by about two orders of magnitude compared to the NR calculation, which gives |ε FS (NR)| = 5 × 10 −13 and 2 × 10 −10 for Ca + and Yb + respectively. This is due to the nonvanishing p 1/2 wave functions at the origin for both Ca + and Yb + as illustrated in Fig. 1 for the latter.
If |ε PS | is smaller than |ε FS |, the experimental search of new force through the IS nonlinearity becomes difficult. This is the case in the shaded regions below dashed lines in Fig. 3. These dashed lines represent the expected sensitivities in experiments of 0.3 Hz precision for Ca + and 1 kHz for Yb + .
In our previous work [7], we described a method to remove the FS nonlinearity by generalizing the linearity relation with three or more transitions. Here we consider an alternative way to avoid the large FS nonlinearity.
The inclusion of two distinct states that have nonvanishing wave functions at the origin, such as s 1/2 and p 1/2 in this work, leads to the large FS nonlinearity as argued in Ref. [7]. This is because the FS is governed by the wave functions inside the nucleus. Hence it is possible to suppress the FS and its nonlinearity by replacing the p 1/2 state with a state whose wave function vanishes at the origin, e.g. the p 3/2 state shown in Fig. 1. Our estimation of the magnitudes of FS nonlinearity for the cases of p 3/2 is presented in the last column of Table 1. We observe that the use of p 3/2 instead of p 1/2 reduces the FS nonlinearity to the level of the NR calculation. This is plausible because of the similarity of the p 3/2 and NR wave functions as seen in Fig. 1. We note that the PS nonlinearity in the mediator mass range below O(0.1) MeV is not significantly affected by substituting p 3/2 for p 1/2 . It is possible in principle to employ a state other than p 3/2 as far as it vanishes at the origin.

Conclusion
We have examined the implication of relativistic calculation in the search of new intraatomic force using the IS nonlinearity. The relativistic wave functions in the singleelectron approximation are obtained by numerically solving the Dirac equation. We have employed the effective potential described by the Thomas-Fermi model with the Helm nuclear charge distribution. We have evaluated the PS and FS nonlinearities with these wave functions. The calculated nonlinearities and the current experimental bounds by the Ca + and the Yb + IS measurements are presented in Table 1 as well as the expectation in future experiments.
In Fig. 2, the bounds on the mass and couplings of the force mediator in the relativistic calculation are compared with those in the NR calculation. We have found that the NR calculation underestimates the magnitude of the PS nonlinearity in the mass range greater than ∼ 1 MeV for Yb + , while the relativistic effect in the PS is less sizable for Ca + .
The current bound and the future sensitivity of IS nonlinearity are also compared with constraints from other experiments and observations in Fig. 3. It is indicated that, although the current bound is weaker than the best other constraints, the IS nonlinearity has a potential sensitivity to probe the unexplored parameter space in future. Future experiments of 1 Hz precision have been supposed in our numerical calculation as an illustration. As mentioned above, IS measurements with the precision of order 1 Hz are likely in the near future. When such data of two or more transitions become available, the sensitivity of IS nonlinearity to new physics is expected to be notably improved.
Moreover, it turns out that the FS nonlinearity is significantly enhanced by the relativistic effect and the new force search beyond the Yb + IS precision of about 1 kHz is limited. The Ca + case is limited at and beyond about 0.3 Hz, so that the problem is less serious.
The large FS nonlinearity is due to the nonvanishing value of the p 1/2 wave function at the origin, unlike the NR p state. Therefore, instead of p 1/2 , we have considered the use of the p 3/2 state, whose wave function vanishes at the origin like the NR p state. We have found that the FS nonlinearity is suppressed in the p 3/2 case as shown in Table 1.
In conclusion, the experimental search of new physics with the IS nonlinearity has a potential sensitivity beyond the existing terrestrial constraints. A proper selection of transitions is important as well as improvements of experimental precision.