Probing new intra-atomic force with isotope shifts

In the development of atomic clocks, some atomic transition frequencies are measured with remarkable precision. These measured spectra may include effects of a new force mediated by a weakly interacting boson. Such effects might be distilled out from possible violation of a linear relation in isotope shifts between two transitions, as known as King's linearity, with relatively suppressed theoretical uncertainties. We discuss the experimental sensitivity to a new force in the test of the linearity as well as the linearity violation owing to higher order effects within the Standard Model. The sensitivity to new physics is limited by such effects. We have found that for Yb$^+$, the higher order effect is in the reach of future experiments. The sensitivity to a heavy mediator is also discussed. It is analytically clarified that the sensitivity becomes weaker than that in the literature. Our numerical results of the sensitivity are compared with other weak force search experiments.


Introduction
One of the main targets of the intensity frontier in particle physics is a new force carrier which is much lighter than the weak scale and very weakly interacts with the Standard Model particles. Such a new particle is examined by low energy experiments and stellar observations.
Among various experiments, measurements of transition frequencies in the development of atomic clocks achieve remarkable precision. Their relative error is expected to be the order of 10 −18 for Yb ion in the near future [1]. Therefore, precision atomic spectroscopy can be considered as a sensitive probe to a new force other than the Coulomb interaction. In contrast to the extreme precision of the experiments, the theoretical calculation of atomic spectra suffers from uncertainties of the many-body system.
A possible way to reduce theoretical uncertainties is to utilize the isotope shift. The atomic spectra differ from one isotope to another. The shifts are considerably small, typically O(GHz), compared with the major part of the Coulomb interaction. The remaining shifts are ascribed to two origins. One is the mass shift and the other is the field shift. The mass shift is caused by the modification of the kinetic term. The field shift depends on the modification of the potential which originates from the nuclear charge distribution differences of isotopes. At the leading order of these effects, the isotope shifts of two different transitions satisfy a linear relation. This linear relation is called King's linearity [2].
The electron-neutron interaction given by a new force carrier also contributes the isotope shift. In general, this effect violates King's linearity. Using this property, the possibility to detect the contributions of the Higgs and Z bosons is studied in Ref. [3]. As revisited in the next section, their effects are highly suppressed by the large scale differences between the electroweak and the atomic physics. However, the violation of King's linearity could be still a good probe of a new force mediated by a boson lighter than about 1 MeV [4,5,6].
In this paper, we study the violation of King's linearity by a new force which is mediated by a light new boson. In particular, its sensitivity to the new physics is discussed by not only numerical but also analytic means. We also study the next-leading-order contribution of the field shift, which also violates the linearity [7,8]. We argue how the field-shift non-linearity limits the experimental sensitivity to the new physics. The expected bounds on the new force in future experiments are compared with the other known constraints. It turns out that King's linearity violation gives us a complementary constraint to a new force in the intra-atomic (sub-keV to sub-MeV) range.
The rest of the paper is organized as follows. In Sect. 2, we revisit the breaking of King's linearity by the higher order effect of the field shift and by the new particle. Afterwards, we present our numerical results of the violations of the linearity in terms of the mass and the coupling of the new boson in Sect. 3. We also compare the future sensitivity of King's linearity violation with other bounds. Section 4 is devoted to our conclusion.

Non-linearity of the isotope shift
We review King's linearity and discuss its violation in this section. Two sources of the linearity violation are examined. One is the next leading order (NLO) of the field shift, and the other is the light new mediator. We call the isotope shift by the exchange of a new particle as the particle shift.
The isotope shift of a transition, which is represented by δ ν, is described by the leading orders of the mass shift and the field shift, and the other contributions as δ ν = Gδ µ + Fδ r 2 + X. (2.1) The coefficients G and F represent the transition-dependent parts of the mass and the field shift, respectively. The last term X stands for the other contributions including the higher-order effects of the above shifts. The reduced mass difference of the isotopes is denoted by δ µ. The transition-independent part of the field shift is given by the difference of the mean square charge radii δ r 2 of the nuclei. Since the precise charge distribution is not clear, to be compared to the reduced masses, we eliminate δ r 2 with another transition. Using the subscripts 1 and 2 for the different transitions, we obtain where We call ω i the modified isotope shift. If the last term is independent of the mass numbers, the modified isotope shifts satisfy the linear relation, that is, King's linearity [2]. In the rest of this section, we firstly study the higher-order correction of the field shift. Since it contributes X, the sensitivity to the new particle is limited by the size of its contribution. Secondly, we discuss the particle shift. The formulation is similar to the field shift. Afterwards, the violation of the linearity is discussed in detail for both of the shifts.

The field shift
The nuclear charge distributions of the isotopes are slightly different from each other. Their small modifications of the potentials are observed as the field shift in the isotope shift.
In the single electron approximation, assuming the spherical symmetry of the system, the field shift of a transition is written as where R N (r) is the radial wave function of the state specified by a set of the quantum numbers represented as N.
The modification of the potential δV (r) is defined as with the charge distribution difference δ ρ(r) between two nuclei of charge Z. The spherical symmetry of δ ρ is the origin of the symmetry of δV . The charge distribution ρ is normalized as 4π dr r 2 ρ(r) = 1. (2.8) Since the total charges of the isotopes are the same, the difference of the potential δV can be expressed as In terms of the perturbation on the charge distribution difference, we do not have to take care the isotope dependence of the wave function at the leading order. The charge distribution difference appears near the origin, and the potential difference does too. The wave function close to the origin is modified from that given by the point charge, which is a good approximation away from the nucleus. Therefore, we separately treat the wave functions in the two regions in calculating the field and the particle shifts. The wave function near the nucleus is evaluated by the power series expansion at the origin, and then it is connected to the solution given by the point charge. The connection point r c is also determined by the smooth connection condition. A more explicit description of our treatment is shown in Appendix A.
For simplicity, we discuss the contribution of a state to the field shift. We write the squared wave function inside the nucleus in a series expansion as where l is the angular momentum. The squared wave function outside the nucleus is denoted by Ξ l . The energy shift is given by the following integration: With Eq. (2.10) and the exchange of the integration order, the first term gives us the Seltzer moment expansion [9], where δ r a = 4π dr r a+2 δ ρ(r). (2.15) The explicit form of r a with the Helm form factor [10] is shown in Appendix B. For l = 0, the leading term of the above expression is proportional to δ r 2 . This term gives the field shift in Eq. (2.1). If the expansion of the wave function can be extended to the outside of the nucleus, the second term of Eq. (2.13) vanishes.
For the case of the Helm distribution, the ratio of the leading and the next leading-order field shift for l = 0 is roughly estimated using the formulas shown in Appendices A  We have found the last approximate formula using r N = 0.519 + 1.00Z 1/3 + 0.103Z 2/3 fm, which is derived by fitting r N in Eq. (B.2) as a function of Z with the mass number A being the standard atomic weight [11]. We have used the atoms of Z ≥ 10 to find the above fit results. For Ca + and Yb + , the NLO-to-LO ratios are 8.69 × 10 −4 and 4.15 × 10 −3 , respectively. The estimate is consistent with Ref. [7].

The particle shift
Since the particle shift is sensitive to only the interaction of electron with neutron, we consider the following potential: where g n and g e are the coupling with neutron and electron, respectively. The new force carrier is supposed to possess the spin s and the mass m. Then, in the single electron approximation, the particle shift is This is just given by the replacement of δV with the Yukawa potential in the field shift. Following the discussion of the field shift, we consider the contribution to the particle shift by a state with the angular momentum l. Omitting the couplings, it is written as where we use the incomplete gamma function, Γ (n, x) = ∞ x dt t n−1 e −t . For a given order of the Seltzer moment in Eq. (2.14), each term of the first summation in Eq. (2.22) is simultaneously eliminated with the corresponding field shift.
If the mass of the mediator is much larger than the inverse of the nuclear size, the particle shift is dominated by the leading term in the first summation of Eq. (2.22). Assuming that the transitions include s-states, the field shift coefficient F is proportional to ξ 0 0 as shown in Eq. (2.14). The particle shift involved in X is dominated by the contribution of ξ 0 0 in Eq. (2.22). Thus Fδ r 2 + X is regarded as the product of ξ 0 0 and the transition-independent (but isotope-dependent) quantity. Eliminating the latter with two transitions results in the linear relation in which the information of the particle shift is lost. This is the reason why the violation of King's linearity is insensitive to the Higgs and Z bosons, as discussed in Ref. [3]. For a light mediator, the integration in Eq. (2.22) gives us the leading contribution. If the mass of the mediator is small enough to cover the whole wave function, this contribution becomes independent of the mass.

Non-linearity
So far, we have obtained some additional contributions to the isotope shift. Including the terms discussed in this section, Eq. (2.1) is modified as (2.23) The higher-order field shift is described byF, and its leading contribution is proportional to δ r 4 . The last term including H is the contribution of the particle shift. These terms have appeared as X in Eq. (2.1), then, in general, they violate King's linearity. In Eq. (2.2), the values F 1,2 are proportional to the differences of the squared wave functions at the origin (see Eq. (2.14)). Omitting the electron spin degree of freedom as the numerical analysis in Sect. 3, three or four distinct electronic orbital states are involved in the two transitions. Here we consider the case that one of the orbital states is an s-state and the others are non-s-states. We observe the following features of the non-linearity: • If both of the transitions include the s-state, its contributions are canceled each other.
• If only one of the transitions includes the s-state, its contribution is suppressed.
In the former case, we obtain X 2 /X 1 = F 2 /F 1 as far as the s-state contributions in X 1,2 are considered, so that the s-state contribution to X 21 in Eq. (2.5) vanishes. In the latter case, assuming that the transition 1 includes the s-state, namely, F 1 = 0 and F 2 = 0, then X 21 = X 2 . That is, the dominant contribution due to the s-state is suppressed. In both of the cases, the non-linearity is induced by higher angular momentum states. We shall use the experimental results of these situations in the numerical analysis later. To obtain effects of the non-linearity by the s-state, the transitions need to include at least two distinct s-states.
The non-linearity given by the particle shift follows the same features. The contributions of the s-states disappear unless the relevant transitions include two or more distinct s-states. Then only the states with the higher angler momenta contribute to the non-linearity. The sensitivity to the new force carrier is scaled with an inverse power of its mass if it becomes much heavier than the typical scale of the wave function. Since the leading contribution is not given by the s-state in the above situation, the scaling is 1/m 4 or worse.
Finally, we mention the extension of King's linearity. The higher-order corrections, such as the NLO field shift above, prevent us from improving the experimental bounds of the light mediator search. However, it is possible to eliminate the effect of the NLO field shift and obtain an extended linear relation if we use additional data of isotope shifts with another transition. The general structure of the linearity is discussed in Appendix C.
If we eliminate the higher-order field shift, the corresponding contributions of the particle shift are simultaneously removed. This means that the above scaling in the heavier mediator region becomes worse. Including distinct s-states helps us to improve the sensitivity in a region between the atomic and the nuclear scales.

Numerical analysis
We show the current status and the future prospect of the new weakly interacting light mediator search with the isotope shift.
The particle shift is measured as the deviation from King's linearity with two transitions of an element. Hence, in addition to the experimental precision, the choices of elements and transitions are also important to measure the effect of the new force. In the following numerical analysis, we compare the calcium ion Ca + and the ytterbium ion Yb + . As shown below, some isotope shifts of Ca + are precisely measured at present and the data satisfy King's linearity within the errors. The relative errors in its isotope shifts are better than 10 −4 . Experiments on Yb + give us similar bounds on the particle shift though the errors are about one order of magnitude worse than those of Ca + . We investigate the future prospect of experimental bounds, simply reducing the experimental error to 1 Hz as an illustration.
The precise wave functions of the states in the transitions are quite difficult to obtain for heavier elements because it is a many-body problem. We calculate them as the wave function of an electron in an effective potential   given by other electrons and the nucleus. The effective potential is estimated by the Thomas-Fermi model, which is a semi-classical approximation of the electrons around the nucleus; see e.g. [12]. Our analysis is done in the non-relativistic limit. In this case, the wave function does not discriminate the spin dependence. The electronic states are characterized by a pair of quantum numbers (n, l). This approximation is good for sand p-states as shown in Appendix D and Ref. [13]. The same data set of Ca + isotope shifts is studied in Ref. [4] using different wave functions from ours. Both results reasonably agree in the sufficiently small mass region.

Experimental data of isotope shifts
The experimental data of the isotope shifts used in our analysis are summarized in Table 1. The masses of isotopes for the calculation of the modified isotope shifts are given in Ref. [11].
Both of transitions include only one s-state, 4s for Ca + and 6s for Yb + . As discussed in the previous section, the s-states do not contribute to the non-linearity of King's plot in these cases. Besides, the field shift of the higher Seltzer moments rapidly become small. Therefore, the p-states numerically dominate the field shift non-linearity.
The constraints to the non-linearity are calculated by the usual χ 2 as described in Appendix E. The bound to g n g e depends on its sign. The weaker bound is employed as the one for its absolute value.

Current experimental bounds and future prospects
The current bounds and the future prospects of King's linearity violation are shown in Fig. 3.2. The lines given by the isotope shifts are the same in both panels. The left and the right panels show other experimental constraints on the scalar and the vector mediators, respectively. The integrand of the particle shift in Eq. (2.20) includes the difference of the squared wave functions. If the mediator is massless, both of the states contribute the integral. Since the integrand flips the sign at a point, the integral disappears around the corresponding mass scale. This cancellation explains the peak structures in Fig. 3.2. The position of cancellation depends on the combination of the wave functions. This means that testing the linearity with various atoms and transitions is important not only to check each against the other but also to exclude the cancellation points.
For the current bounds by the isotope shifts, Ca + and Yb + give us similar bounds, although the experimental errors of Yb + is about one order of magnitude worse than those of Ca + . Accordingly, the sensitivity of Yb + is about one order of magnitude better than that of Ca + in the prospected bounds with the error of 1 Hz.
The field-shift non-linearity of Ca + and Yb + appear at 1.1 × 10 −2 Hz and 4.7 Hz, respectively. These frequencies are interpreted in terms of the mass and the coupling of the light mediator as indicated by the dashed lines in Fig. 3.2. Once the experimental sensitivity reaches the line, the bound on the particle shift is not improved without the further subtraction of the field-shift non-linearity as discussed in Appendix C. The Yb + line of the expected sensitivity lies below the field-shift non-linearity indicated by the dashed line, while the Ca + line does not. If the precision of the Ca + measurement is improved, it may cover the region of smaller coupling. However, in this case, the non-linearity of the NLO mass shifts, which are more significant for lighter elements, also ought to be taken into account.
The experimental bounds are obtained as the allowed size of the non-linearity. They depend on the sign of the coupling with the light new mediator. Hence, the constraints to the mediator are changed at the peaks since the sign of the particle shift is changed there. As already mentioned, we have adopted the sign giving the weaker constraint to calculate each of the current experimental bounds. Thus, the current bounds shown in Fig. 3.2 are the conservative ones. Now, we turn to other constraints by light particle searches in the same mass region. We consider the terrestrial bounds which are obtained as the product of the individual bounds on the couplings with electron and neutron. The constraint of the electron coupling is mostly given by the electron g − 2 [17]. The region above 20 MeV is bounded by Babar [18]. The beam dump experiments strongly constrain the region from 100 keV and 1 MeV to about 10 MeV for the scalar and the vector mediators, respectively [19,20]. As summarized in Ref. [5], the experimental bounds on the neutron coupling is given by the combination of several low energy neutron experiments [21,22,23]. The obtained bounds are much stronger than the current bounds by the isotope-shift non-linearity. However, the future prospects of the sensitivities are better than the above terrestrial bound in the region between about 10 eV and 1 MeV.
The mass region less than about 100 eV is covered by the fifth force search [26,27]. Since this bound is very strong, the mediator lighter than about 100 eV is almost excluded for the entire coupling region plotted in Fig. 3.2.
The stellar cooling bounds are also very strong in the region lighter than about 1 MeV [24]. However, if the coupling becomes large, the light new boson cannot take the energy to the outside. It is not sure that the stellar cooling observation can constrain the boson for sufficiently strong couplings. The brown dotted line in Fig. 3.2 indicates this limitation given by Ref. [25]. We consider the isotope-shift non-linearity as well as other terrestrial experiments to be useful since it provides us constraints independent of the stellar dynamics.
The mass and the coupling of the light vector boson suggested by the Atomki 8 Be experiment [28,29] are also shown in the right panel of Fig. 3.2. The result in Ref. [4] is that the region can be excluded by the future Yb + measurement. However, our result indicates that the sensitivity cannot reach the coupling of the vector. Even worse, the field-shift non-linearity prevents us from probing the parameter region. We note that the present work brown regions below about 500 keV is constrained by the stellar cooling bounds given by Ref. [24]. As described in Ref. [25], the stellar bounds have uncertainty above the brown dotted line. The shaded gray regions below about 100 eV are restricted by the fifth force experiment [26,27]. The black line in the right panel stands for the region indicated by the Atomki anomaly [28,29]. and Ref. [4] employ different transitions of Yb + , and our result is derived from the existing data of the Yb + isotope shifts.

Conclusion
The isotope shifts can be precisely measured with the technique developed in the atomic clock experiments. In the measurements, the weakly interacting light force carrier can be tested as the violation of King's linearity. We have investigated two sources of the non-linearity. One is the higher-order effect of the field shift and the other is the above effect of the light mediator, the particle shift. Because of the larger particle shift of heavier elements, the sensitivity to the light mediator of Yb + is about one order of magnitude better than that of Ca + . However, Yb + also possesses the sizable field-shift non-linearity. In the future prospects with the 1 Hz error, the sensitivity of Yb + reaches the non-linearity given by the field shift. Then its sensitivity to the light mediator is not improved without the further development of the formulation of the linearity. On the other hand, the field-shift non-linearity of Ca + is about one order of magnitude smaller than that of Yb + in terms of the coupling of the light mediator.
The future prospect gives us the bounds better than the other low energy experiments in the region from 100 eV to 1 MeV. This region is also bounded by the observation of the stellar cooling. The terrestrial experiments including the isotope-shift non-linearity tell us the complementary information.
In our study, the scaling of the sensitivity is analytically clarified as 1/m 4 in the region where the force carrier is heavy. In the numerical analysis, this behavior is indeed observed if the boson is heavier than about 100 keV. As a result, the expected exclusion lines with the error of 1 Hz cannot reach the parameter region of the Atomki anomaly. Furthermore, the favored region cannot be proved because of the field shift non-linearity even if the experimental precision is improved.
Our analysis uses the non-relativistic limit, the mean field approximation with the Thomas-Fermi potential, and the assumption of the Helm distribution. It is required to improve the theoretical methods for more precise and extensive studies of the isotope-shift non-linearity including other atoms and transitions.

A. Wave function inside the nucleus
The Schrödinger equation of the radial direction is Since we consider a bound state, the potential V (r) and the energy E are negative in this equation. We are interested in the wave function near the origin. The potential is expanded as where v 0 < 0, and v 1 = 0 for a nuclear charge distribution without a cups at the origin such as the Helm distribution described in Appendix B. It is supposed that the wave function is also given as the following series: Then the Schrödinger equation is expressed as The first term leads to χ l 1 = 0. The entity within the large parentheses gives us the recurrence relation to determine the coefficients χ l i . The first equation of the above relation with i = 0 is Using E/v 0 1, we obtain The rest of the coefficients are recursively obtained. We use the wave function of O(r 2 ). The squared wave function inside the nucleus is This wave function is smoothly connected to the wave function with the point charge, i.e., the solution outside the nucleus. The coefficient χ l 0 and the connection point r c are both obtained by the condition to connect. It is convenient to express the wave function with the point charge as r l (a 0 + a 1 r) near the nucleus. Then the above quantities are analytically calculated as Strictly speaking, the wave function inside the nucleus varies from one isotope to another because v 0 depends on r N as shown in Eq. (B.13). Since this is a higher-order effect, we use the fixed r N of A = 44 for all the Ca isotopes and A = 173 for the Yb ones.

B. The nuclear charge density and the potential with the Helm distribution
The Helm distribution of the nuclear charge [10] is defined as In our numerical analysis, we use the following parameters given by Ref. [30]: This distribution is given by smearing of the uniform density distribution as shown below. The limit of s → 0 corresponds to the uniform distribution with the radius of r N . The spatial density is obtained: where the error function Erf(x) is defined as The distribution is also expressed as where θ is the step function and the smearing function ρ g is given by the Gaussian distribution as We can check the whole space integration of this density indeed becomes unity, namely, The electrostatic potential induced by the above density is This function is even with respect to r. The potential at the origin is (B.13) In the Seltzer moment expansion, we need the mean values of r n , r n =4π dr ρ(r)r 2+n (B.14) where we introduce the confluent hyper geometric function of the first kind, The relevant terms in our calculations are

C. Generalization of the linearity
Eliminating the difference of the mean squared radii with two distinct transitions, King's linearity is obtained. It turns out that the higher-order field shift and even the particle shift can be similarly eliminated with additional transitions. We formulate this procedure and derive the generalization of the linearity as follows.
For simplicity, we discuss it with the higher-order field shifts at first. The isotope shift is written as where the terms including F (k) i are the higher-order field shifts in terms of the Seltzer moment. The last term X i stands for the other non-linearity. If we have a sufficient number of transition data, the above expression can be rewritten as For convenience, we write the matrix in the right hand side as T . If T has an inverse, we multiply it to the above equation. The difference of the reduced masses δ µ is precisely measured, while the mean radii are not.
Multiplying T −1 to the above equation, the first element is Dividing this equation by δ µ, we obtain the generalized linear relation which is free of the field-shift non-linearity. If X i = 0 for any i, the data of the modified isotope shifts with n different transitions are on an n − 1 dimensional plane. The above procedure is easily extended to the other non-linearity. Firstly, the given term is separated into the wave-function-dependent part and -independent part. This way of separation is not unique, however, the result is independent of this detail. Secondly, the wave-function-dependent part is embedded in the above matrix T , and then we multiply the inverse. Finally, the precisely measured element, like δ µ above, gives us the linear relation we want. This means that we may simultaneously obtain several linear relations.
For example, we consider the elimination of the non-linearity by the particle shift. In order to do so, the factor independent of the wave function (A − A) is appended to the above vector, and the other part is embedded in T . Then a new linear relation is derived in the same manner as above. Since the factor (A − A) is also determined with a high precision, the additional linear relation is found from the different element. Even if King's linearity violation is observed, these generalized linear relations are preserved as long as the non-linearity originates from the particle shift.

D. The Thomas-Fermi potential
The Thomas-Fermi model is a semi-classical approximation of the electrons around the nucleus; see e.g. [12]. The electrons are regarded as the free fermion gas of zero temperature. Two electrons occupy the phase space of (2πh) 3 from the bottom. The self-consistent potential is represented by the universal function so called the Thomas-Fermi function.
This function is the solution of the following non-linear differential equation: One boundary condition is given at the origin as χ(0) = 1. For positive ions, the other boundary condition is imposed to satisfy  where χ(x) vanishes at x 0 , and n is the positive charge of the ion. Since we consider an electron in a singly charged positive ion, the mean potential is given by the other electrons and the nucleus, namely, n = 2. The resulting potential energy is where x = 4 3 2Z/9π 2 m e αr. The boundary x 0 is approximately given by x 0 = −8.964 + 7.341Z 1/3 . The wavelengths of the relevant transition evaluated in the Thomas-Fermi potential are shown in Table 2 as well as the corresponding experimental data.

E. Statistics
We use the following formulae in the numerical analysis. The data of the modified isotope shifts of two transitions are denoted by x a and y a , and their standard deviations are σ xa and σ ya , respectively. The subscript a indicates an isotope pair. The term violating the linearity is represented by εs a . The parameter ε stands for the wave function independent part, e.g., the coupling of the particle shift. The χ 2 of the fit function can be written as These conditions show us the χ 2 minimum, which is zero in the present work because of the additional parameter ε. We calculate the χ 2 minimum as a function of ε. Writing it as χ 2 ε , the fit quality is measured as χ 2 ε /dof, where dof = 1 in our analysis.
The extension to the case of more than three isotope pairs or more than two transitions is straightforward.