Probing new light force-mediators by isotope shift spectroscopy

In this Letter we explore the potential of probing new light force-carriers, with spin-independent couplings to the electron and the neutron, using precision isotope shift spectroscopy. We develop a formalism to interpret linear King plots as bounds on new physics with minimal theory inputs. We focus only on bounding the new physics contributions that can be calculated independently of the Standard Model nuclear effects. We apply our method to existing Ca+ data and project its sensitivity to possibly existing new bosons using narrow transitions in other atoms and ions (specifically, Sr and Yb). Future measurements are expected to improve the relative precision by five orders of magnitude, and can potentially lead to an unprecedented sensitivity for bosons within the 10 keV to 10 MeV mass range.


I. INTRODUCTION
The Standard Model of particle physics (SM) successfully describes multiple observations up to the TeV scale, and is theoretically consistent up to a much higher energy. However, the SM cannot be a complete description of Nature. For example, it lacks a viable dark matter candidate and can neither explain the observed matterantimatter asymmetry of our Universe nor neutrino oscillations. In addition, the SM suffers from hierarchy issues both in the Higgs sector and the fermionic sector. These experimental observations require new physics (NP) beyond the SM, however, none of these observations point towards a specific new theory or energy scale.
The quest for NP is pursued in multiple directions. Current efforts with colliders such as the LHC form the energy frontier, probing directly the TeV energy scale. Other accelerators, such as B-factories, NA62 and neutrino experiments, form the intensity frontier that broadly probes the MeV-GeV scale. Atomic physics tabletop experiments form a third frontier of precision measurements (see e.g.: [1][2][3][4][5], for a review see [6][7][8]) where sub-MeV physics can be efficiently tested. It is interesting to note that NP that may account for the hierarchy issues could be new light scalars that couple to matter fields [9][10][11][12][13][14][15]. To convert the high precision offered by atomic and molecular spectroscopy into sensitivity to fundamental new physics, one either has to acquire similar theoretical accuracy of atomic structure or alternatively seek for unique observables that are insensitive to theoretical uncertainties.
In this paper we show that precision isotope shift (IS) spectroscopy may probe spin-independent couplings of light boson fields to electrons and neutrons. The idea is to extract constraints from bounds on nonlinearities in a King plot comparison [16] of isotope shifts of two narrow transitions [17]. We develop a new formalism to interpret these measurements in the context of searching for new light force carriers and propose several elements and transitions that can be used for such analyses. We recast existing measurements into bounds and provide an estimation for the sensitivity of future measurements, see Fig. 1. The validity of our method to bound NP does not rely on the knowledge of the SM contributions to King plot nonlinearites. Its constraining power, however, is limited by the size of the observed nonlinearities. In case that Kings linearity is established, at the current state- σi specified in the labels). Constraint from existing IS data: Ca + (397 nm vs. 866 nm [18], solid red line). IS projections (dashed lines) for Ca + (S → D transitions), Sr + , Sr/Sr + , and Yb + . For comparison, existing constraints from other experiments (shaded areas): fifth force [19,20] (dark orange), (g − 2)e [21,22] combined with neutron scattering [23][24][25][26] (light blue) or SN1987A [27] (light orange), and from star cooling in globular clusters [28][29][30] (orange). The gray line at 17 MeV indicates the yeyn values required to accommodate the Be anomaly [31,32].
of-the-art experimental precision, and baring cancellation between the SM and NP contributions, world-record sensitivity in a certain mass range will be achieved.

II. FACTORIZATION OF NUCLEAR AND ATOMIC EFFECTS IN ISOTOPE SHIFTS
We now discuss the scaling and factorization properties of IS which we use to probe NP in this work. Consider an atomic transition, denoted by i, between narrow atomic states. The difference in the transition frequency ν comparing the isotopes A and A is the IS, At leading order (LO) the IS receives contributions from two sources, mass shift (MS) and field shift (FS). Mass shift arises due to a correction to the kinetic energy of atomic electrons due to the motion of the nucleus. For independent electrons, this is just replacing m e by the reduced mass but if electrons are correlated, this could be orders of magnitude larger. Field shift originates from different contact interactions between electrons and nuclei in isotopes. Putting these two leading contributions together, the IS can be phenomenologically written as where the two terms represent MS and FS respectively [16,33]. We define µ AA ≡ m −1 A − m −1 A , where m A and m A are the masses of isotopes A and A .
The quantity δ r 2 AA is dominated by the difference in the mean squared charge radii of the two nuclei but can include other contact interactions. Both µ AA and δ r 2 AA are purely nuclear quantities that do not depend on the electronic transition i. Note, however, that µ AA is known with high precision, whereas δ r 2 AA is known only to a limited accuracy. The parameters K i and F i are isotope-independent, transition-dependent coefficients of the MS and FS, and their precise values are unnecessary in the observable we construct. Each term of Eq. (2) is a product of a purely nuclear quantity and a purely electronic quantity, resulting in the factorization of nuclear and electronic dependence. This is known as LO factorization.
Given two electronic transitions, i = 1, 2, one can eliminate the uncertain δ r 2 AA giving a relation between the isotope shifts ν AA 1 and ν AA 2 . In terms of the modified with Equation (3) reveals a linear relation between mν 1 and mν 2 , giving rise to a straight line in the so-called King plot of mν 2 vs mν 1 [16]. It is important to stress that the linearity of this equation holds regardless of the precise values of the K i and F i electronic parameters. Testing linearity necessitates at least three independent isotope pairs in two transitions, which constitutes a purely data driven test of LO factorization.
The formulae in our treatment of NP will be simplified greatly by introducing a geometrical description of LO factorization. As we will now explain, King linearity is equivalent to coplanarity of vectors. For each transition i, we can form a vector The nuclear parameters of field and mass shift, µ AA and δ r 2 AA can also be written as vectors −→ mµ and − −−− → mδ r 2 in the same space (notice that −→ mµ ≡ (1, 1, 1)) and hence Eq. (2) becomes 1 Below we will adopt the notation of adding an m to "modified" (i.e. normalized by µ AA ) quantities, such as mδ r 2 AA ≡ δ r 2 AA /µ AA .
In this language LO factorization implies the following qualitative statement: any vector of reduced isotope shifts, − → mν i , must lie in the plane that is defined by −→ mµ and − −−− → mδ r 2 , as illustrated in the cartoon in the left panel of Fig. S1.
Note that, because the direction of − −−− → mδ r 2 in this space is uncertain, theory does not tell us in which direction this plane is oriented. However, by measuring two IS vectors, − → mν 1 and − → mν 2 , we can test this statement by asking whether the three vectors − → mν 1 , − → mν 2 , and −→ mµ are coplanar. The coplanarity of these vectors corresponds to King linearity as we can see by rewriting Eq. (3) in vectorial form − → mν 2 = K 21 −→ mµ+F 21 − → mν 1 . Like King linearity, coplanarity is a purely data driven test of LO factorization since it is independent of theoretical input. A change in K i and F i will merely change the direction of − → mν 1 and − → mν 2 within the plane, but the qualitative statement of coplanarity remains.
In this vector language we can provide a compact expression for a nonlinearity measure, In terms of the King plot, NL is the area of the triangle spanned by the three points shown in Fig. S2. Equivalently, in our geometrical picture it is the volume of the parallelepiped defined by − → mν 1,2 and −→ mµ. A given data set is considered linear if NL is smaller than its first-order propagated error σ NL = Σ k (∂NL/∂O k ) 2 σ 2 k where the sum runs over all measured observables O k (modified frequency shifts and isotope masses) with standard deviations σ k .

III. NEW PHYSICS AND VIOLATION OF KING LINEARITY
We now include a NP contribution by adding a third, also factorized, term to Eq. (2), namely X i depends on the form of the new potential and on the electronic transition, while γ AA depends only on the nuclear properties. The parameter α NP is the NP coupling constant which we would like to probe. Let us first mention two cases of NP which we do not expect to be able to probe by testing King linearity. For short-range NP (shorter than the nuclear size), the electronic parameters X i will be proportional to those of FS, X i ∝ F i . In this case the NP term can be absorbed by redefining δ r 2 AA . Also, if the new physics couples to electrons and nuclei according to their electric charge (such as the case of dark-photon [34]), γ AA = 0. There may also be cases in which NP can accidentally be absorbed by redefining F i . However, a long-range force with couplings not proportional to the electric charge (and barring an accidental cancellation) can be severely constrained by tests of King linearity.
Equation (3) written in vectorial form becomes where h is the NP vector in reduced frequency units, that is h AA ≡ γ AA /µ AA and X 21 ≡ X 2 /X 1 . One can see that NP can lead to a deviation from coplanarity if and only if (i) the new force is not short-range, X 21 = F 21 ; (ii) h is not aligned with any linear combination of −→ mµ, − → mν 1 or − → mν 2 . By solving the set of equations (7) one finds an expression for α NP that is needed to yield a particular dataset assuming NP is the dominant contribution to nonlinearity. If linearity holds then α NP σ αNP = Σ k (∂α NP /∂O k ) 2 σ 2 k . Hence, the sensitivity to probe α NP is lost in the limit where the denominator in Eq. (9) vanishes, because the NP contribution to nonlinearity is It is straightforward to check that this happens under the conditions specified below Eq. (8).
The presented method of limiting α NP , Eq. (9), contains theory input only in X i and h AA which describe how NP affects the IS. The SM contribution in the factorized limit is fully parametrized by the observables ν i and µ. The form of h AA depends on the assumed couplings of new physics to nuclei. For example, if the new interaction couples to quarks, then we expect that h AA ∝ AA [17,35]. The atomic transition-dependent factors X 1,2 can be reasonably calculated by a manybody simulation (see the next section). This strategy is analogous to a search for NP, say, at the LHC, where all SM backgrounds are estimated using data driven methods and Monte Carlo simulation is used only in estimating the signal cross section.
Thus far, most measurements of scalar-isotope shifts have been consistent with King linearity (see, however, the case of Samarium [36]). Nevertheless, some level of nonlinearity is expected to arise from SM higher-order contributions [37][38][39][40]. These contributions, that are related to nuclear physics and electronic-structure dynamics linked together, are presently not understood in a quantitative manner for many-electron systems. One possible source of nonlinearities is of the form of a field shift that depends on the isotope mass. Precision calculations recently showed that this effect is of O(10 −3 −10 −4 ) in light atoms [41]. Likewise, such contributions in heavier elements with Z = 20 − 87 [39], but only for S → P transitions, are estimated to be of a similar order. Hence, matching the precision of future measurements motivates the calculation of the remaining higher-order corrections.
If a deviation from King linearity is observed, it will be difficult to distinguish the NP and SM contributions to the nonlinearity. In this case there are two options in which further insight on NP can be obtained. The first requires that the theory of King nonlinearity would advance and enable us to subtract the SM contributions, and in the process possibly gain new insight on the nature of nuclear effects in IS. To add to that, since nonlinearity in the case of NP is universal and in the case of SM specific to particular atomic configurations, a comparison between measurements in different systems will be beneficial. The second relies on the fact that NP forces are of longer range than nuclear effects which require overlap of the electronic wavefunction with the nucleus. Hence it might be possible to identify an observable that is less affected by the nucleus, but is still sensitive to the presence of long-range new physics interactions. In this regard, IS measurements involving Rydberg states might provide a smoking gun for the above types of NP.
For the proposed method to be effective, the element and the specific transitions should be chosen carefully. First, to make a significant progress as compared to current precision, we consider narrow optical clock transitions. The most accurate frequency measurements to date, with a relative error of 10 −18 corresponding to sub-Hz accuracy, have been performed on narrow opticalclock transitions in laser-cooled atoms or ions [42][43][44][45][46][47]. Second, since the hyperfine interaction of electrons with the nucleus is a source for King nonlinearity [37], we consider only even isotopes without nuclear spin.

IV. CONTRIBUTION OF NEW BOSONS TO ISOTOPE SHIFTS
In this section we discuss how theoretical IS predictions are modified in the presence of hypothetical new force carriers of spin s = 0, 1 or 2 and mass m φ which couple to electrons and neutrons with strength y e and y n , respectively. The effective spin-independent potential mediated by such bosons between the nucleus and its bound electrons is V φ (r) = −α NP (A − Z)e −m φ r /r, where α NP = (−1) s y e y n /4π. Note that NP could also couple to protons, though without affecting the linearity of the King plot, hence we neglect such a coupling here.
To calculate the effect of this NP potential on atomic energies we use the "finite field" method where the potential is added directly to the Dirac equation in our many-body computations. The atomic structure calculations are variants of the combination of configuration interaction and many-body perturbation theory (CI+MBPT) [48]. For the single-valence electron ions Ca + and Sr + , we create an operatorΣ (see for example [49]) representing core-valence correlations to second order in the residual Coulomb interaction. This operator is added to the Dirac-Fock operator, along with the NP potential, to generate self-consistent solutions. In this approach, the sensitivity of a transition i between electronic states a and b (i = a → b) can be expressed where ab is the difference of the energy levels of the states a, b, evaluated as a function of α NP and the derivative is taken numerically at α NP = 0. For neutral Sr, which has two valence electrons above closed shells, we use the CI+MBPT method as described in [50]. Briefly, we find the self-consistent solution of the Dirac-Fock equations, including the NP potential, for the closed-shell core (i.e. the V Ne−2 potential where N e is the total number of electrons). In this potential we generate a set of B splines [51,52] which form a complete basis set. Valence-valence correlations are included to all orders using CI, while the core-valence correlations are included using second-order MBPT to modify the radial integrals. The Yb + case is more complicated because of the hole transition, 4f 14 6s → 4f 13 6s 2 . For this ion we use the particle-hole CI+MBPT method in the V Ne−1 potential [53] which has previously been used for Hg + .
The many-body calculations can be cross-checked by perturbation theory, which yields where |Ψ(r)| 2 is the electron-density evaluated in the absence of NP, and h AA = AA amu for the NP contribution in Eq. (7). As a cross-check of our many-body calculation, we use GRASP2K [54] to evaluate |Ψ(r)| 2 and compute X i using Eq. (12) for several Ca + transitions. We find good agreement between the two methods. We identify three regions of NP interaction range, separated by the electron wavefunction size, a 0 /(1+n e ), and the nuclear charge radius, Here a 0 ≈ (4 keV) −1 is the Bohr radius and n e is the ionization number. For m φ (1 + n e )/a 0 , the "massless limit", the interaction range is larger than the atomic size and V φ ∝ 1/r so that X i becomes independent of m φ . For intermediate masses, (1 + n e )/a 0 m φ 1/r N , the interaction range is within the size of the electron wavefunction, and the potential V φ ∝ e −m φ r /r is massdependent. Hence, detailed knowledge of the electronic wavefunctions is necessary to evaluate the effect of NP. In the heavy mass limit, m φ 1/r N , the interaction range is shorter than the nuclear radius and V φ ∝ δ(r)/(m 2 φ r 2 ). In this limit, the NP and nuclear charge-radius effects are approximately aligned since This results in a suppressed sensitivity for new physics which scales as (X 21 − F 21 ) → 0, see Eq. (3), and [17].
In the massless limit, X i can be estimated without a detailed computation of the atomic wavefunctions, as in this case the effective potential is Coulomb-like and thus its effects are approximately accounted for by a shift of α, see Appendix II. We do not estimate the bounds on α NP in the heavy mass limit as in this limit NP effects are indistinguishable from those of finite nuclear size. Bounds are therefore suppressed by a factor of O(r N /a 0 ).

V. CURRENT BOUNDS AND PROJECTIONS
Here we derive the constraints on the product of electron and neutron coupling, y e y n , from existing IS data of Ca + and project the bounds for different transitions alkali-like systems in the 10 eV-50 MeV mass range, assuming that better IS data will be available in the future. Our results are summarized in Fig. 1 and Tab. I.

A. Constraints from King linearity
We apply our method to available IS data of Ca + (solid line of Fig. 1). In the massless-boson limit, m φ 10 keV, the bound is essentially independent of m φ . At the high mass limit, we expect that F 21 = X 21 . Since the theoretical control of F 21 is worse than the experimental error, one can get an incorrect m φ dependence of the y n y e bound at that limit. However, the ratio F th 21 /X 21 (F th 21 is the theoretical value calculated in the absence of NP) has much smaller error. Thus, in order to account for the reduction in sensitivity as m φ increases, we rescale the y e y n bound by ( is the measured value of F 21 . We verified with GRASP2K that this factor does not change by more than a few percent if the charge radius is changed by order one (that is known to a few percent accuracy, hence this is a rather conservative approach). Indeed we see that for m φ > Zαm e the limits get weak, and the sensitivity decreases approximately as m −3 φ for large masses. In Appendix IV we give two heuristic arguments that obtain this asymptotic scaling of our loss of sensitivity 2 : the first is based on approximating Eq. (12); and the second is based on a non-relativistic QED (NRQED) effective theory approach.
Among the non-IS experiments that probe y e y n − m φ parameter space, we consider here only the ones most sensitive to new light bosons coupled to electrons and neutrons. The shaded regions in Fig. 1 summarize the current reach of these experiments. We stress, however,  that some of them are derived involving of further theory assumptions, in contrast to our method which relies on few theory inputs. For new bosons lighter than few×100 eV, fifth force experiments [19,20] are potentially sensitive. Since the interaction range covered by these experiments is much larger than the atomic size, only forces with non-zero atomic coupling can be probed. For illustration we show in Fig. 1 the fifth force bound applicable to U(1) B−L gauge bosons [56]. Furthermore, separately y n is constrained by various neutron scattering experiments [23][24][25] and y e by the anomalous magnetic moment of the electron (g − 2) e [21,22] and by electron beam-dump experiments for m φ > 1 MeV.
Both y e and y n are also severely constrained by globular cluster energy loss for masses m φ 10 keV [27][28][29][30] down to y e y n < 10 −25 and y e by sun cooling [57,58]. Couplings to nucleons in the 10 −10 − 10 −7 range for m φ 100 MeV may be also excluded by energy loss in the core of SN1987A [27,59]. However, such astrophysical bounds might be avoided in certain models such as chameleon, see [60] and references therein. In order to derive an upper bound on y e y n , we combine for each mass the best constraint on y n from neutron experiments with y e either from (g − 2) e or from astrophysics.

B. Prospect for future measurements
As the precision of optical spectroscopy continues to improve, higher accuracy IS measurements in different systems can be achieved in the near future. Accordingly, we estimate the sensitivity that would be achieved for several transitions in alkali or alkali-earth ions or atoms, given the improved accuracy.
Here we consider a comparison between the two finestructure split electric quadrupole transitions in Ca + and Sr + . A comparison between the optical clock transitions in Sr + and Sr, and the quadrupole and octupole transitions in Yb + are also presented. In principle, to enhance the sensitivity of our method, it is desired to compare transitions that involve levels that are as different as possible. For this reason comparing the two fine-structure split electric quadrupole transitions in Ca + or Sr + is not ideal, especially when compared to the sensitivity of the E2 and E3 lines in Yb + or comparing the E2 line in Sr + with the intercombination line in Sr. We include these transitions in our projections since their high-resolution IS measurement is experimentally simpler.
All the transitions above are expected to be measured with 1 Hz accuracy. Under the assumption that King linearity will hold in those future measurements and following Appendix III, the projected bounds are plotted in Fig. 1 (dashed lines) as a function of m φ and summarized in the lower part of Tab. I. The resonance structures, around the 10 keV scale, arise from cancellations in the denominator of Eq. (9). These local losses of sensitivity at different masses per atomic system provide another motivation for IS measurements in complementary systems for a good coverage of the parameter space.
The various projections with 1 Hz accuracy significantly improve the bounds in the m φ ≥ 10 keV region in parameter space. For lower m φ they are weaker than astrophysical bounds. However, astrophysical bounds are subject to large uncertainties and can be broken by models such as the chameleon effect [60]. Thus, an independent laboratory bound in this low-mass region is nevertheless worthwhile. For m φ ∼ a few MeV the projections of Ca + (S → D transitions) and Sr + are comparable to y e y n from neutron scattering [24] and (g − 2) e . Since neutron experiments are affected by uncertainties [25,[61][62][63] such as those in electron-neutron scattering length, nuclear input values and missing higher-order terms in the neutron-scattering cross section, the bounds in the high-mass range well above the neutron energies of E n < 10 keV [24] should be understood as an indication of the order of magnitude. Consequently, theoretically cleaner IS probes at the same order will already improve the bound robustness. Note that a Sr/Sr + (Yb + ) IS comparison would become more effective than other existing methods in probing new bosons above ∼ 10 keV already with 100 Hz (1 kHz) accuracy (the bound related to Sr/Sr + constructed from comparison of transition involving neutral and ion systems suffers from some numerical instabilities for masses above 20 MeV and thus is not shown). Finally, let us note the range of y e y n needed to explain the Be anomaly [31,32] can be probed by future IS measurements of Yb + at the 1 Hz level.
Here we estimate the NP contribution X i to IS in the special case where the force-carrier is much lighter than the inverse atomic size, m φ (1 + n e )/a 0 ∼ O(few keV). Since in this limit the effective potential is Coulomb-like, V φ (r) (A−Z)α NP /r, X i can be simply estimated through a shift of the fine-structure constant α, without a detailed calculation of the electronic wavefunctions. At fundamental level, the Coulomb potential is modified by the shift In the absence of NP, the binding energy of the atomic level a for isotope A scales as where Z a eff ≡ Z − σ a is the effective nuclear charge seen by the valence electron in the state a, and I A a is a constant independent of the charge (modulo O(αZ a eff ) 4 corrections from the fine-structure). The constant σ a > 0 accounts for the screening due to inner electrons. A similar screening effect may occur for the new physics force such that Eq. S4 implies to shift the physical observables as where the constant σ a accounts for the screening of the nuclear NP charge by inner electrons. Note that precise knowledge of this constant is not crucial since it is universal (as a first approximation) for all isotopes and will therefore cancel in X i . Hence the prediction for the IS of the transition i = a → b is (in natural units, i.e. the reduced Planck constant is set to = 1) and expanding to leading order in α NP α yields Since IS are typically orders of magnitude smaller than the transition frequencies (with the exceptions of very degenerate states such as in dysprosium [64,65]), we can take E A ≈ E A in the NP contribution above and matching to Eq. (7) we find γ AA = A − A and By a combination with Eq. (9), the above expression provides a reasonable estimate of the constraint on α NP without the need for an accurate knowledge of electronic densities. For example, we found that the upper bounds on y e y n obtained by evaluating X i with Eq. (S9) and with the CI+MBPT method are comparable and only differ by O(1) factors.

III. PROJECTING FUTURE BOUNDS
Our procedure above applies to cases with enough experimental data. For systems lacking (sufficiently precise) measurements, we can still derive projections provided that an acceptable estimation of the F 21 constant is available from either theory calculation or hyperfine splitting data (whenever available). Assuming the observation of linearity and global experimental uncertainties σ i of ν AA i , the only missing information is how much the new physics vector, h, points towards the linearity plane. While a precise determination of the vector component requires data, the projection of the NP vector h AA = AA amu along the mass shift direction is easily obtained without the need of experimental input. Therefore, a best-case projection 3 [σ αNP ] proj can be obtained by neglecting the possible additional alignment of NP with nuclear effects, where ∆A min(max) j ≡ min (max)[A − A j ] and σ i is the assumed standard deviation of IS measurements in transition i. Note that the sensitivity to α NP is weaker by a factor of A/∆A min j than the naive expectation of Refs. [17,66,67]. The reason is that the NP physics vector lies mostly in the linearity plane, in particular it has a large projection along the mass shift direction. Eq. (S10) implies that elements with small A/(∆A min j ∆A max j ) are preferred. However, if the mass shift dominates over the field shift, which is the case for light elements, the sensitivity is reduced as well.

IV. SCALING OF HIGH MASS LIMITS
As we discussed in Section III, King linearity (or coplanarity) is only sensitive as a probe of long-range forces, which extend beyond the nuclear size. This is because the contribution of a short-range interaction to the isotope shift can be absorbed into the contact interactions such as the nuclear charge radius. From the perspective of an atom, a new interaction begins to approach a contact interaction when its range is shorter than the effective radius of the inner most K-shell electron, (Zαm e ) −1 = a 0 /Z. We thus expect our method to start losing power for mediators heavier than Zαm e as can be seen in Fig. 1.
It is instructive to investigate like what power of m φ /(Zαm e ) one would expect the limits to decrease. As shown in Eq. (8) and the surrounding discussion, the contribution of NP to nonlinearity is proportional to (X 2 − X 1 F 21 ). When the NP contribution is aligned with the FS, X i ∝ F i , this factor vanishes. We must thus investigate how X i and F i differ for the various transitions in the limit of m φ → ∞. In our proposed procedure, F i is measured from data and X i is estimated from theory, as described in Section IV, leading to the limits shown in Fig. 1. Since these two quantities are extracted using different methods we would like to ensure that the alignment of X i and F i is captured correctly and that the weakening of the bounds for m φ > Zαm e agrees with our theoretical expectation.
The asymptotic behavior of the limits in Fig. 1 exhibits an approximate m 3 φ behavior. We will now show that this may be expected on theoretical grounds, first in position space, using a hydrogen-like approximation and then in momentum space, using an effective field theory (EFT).
To leading order in the small new physics coupling, the contribution of a new Yukawa potential to the energy of the atomic level a is To study the scaling behavior in a simple case, we approximate the multi-electron crudely as a single-electron atom with an effective Z (a) eff which will account for the screening of the nucleus by the inner shells. Of course, the full calculation of Section IV accounts for the multi-electron effects fully and the approximation we use here should be taken as a toy model to study the scaling of our bounds. We consider atoms in the nS state, with (n = 1, 2, . . .), which is relevant for us since many of our isotope shifts involve an S state. Using the explicit form of non-relativistic hydrogen-like wavefunctions, it is simple to show that in this case where we assumed that m φ Zαm e and expanded in powers of 1/m φ , keeping the first two terms. The leading term is to be expected, and is identical to the shift in energy from a contact interaction potential δ 3 (r)/m 2 φ . This term will thus lead to isotope shifts in the a → b transition that are proportional to |Ψ a (0)| 2 − |Ψ b (0)| 2 and thus proportional to the FS. This term therefore cancels in the combination (X 2 − X 1 F 21 ). The second term in Eq. (S12) has an additional factor of Z eff leading to an IS proportional to Z  We can also understand the scaling of our bounds at high mediator mass in an EFT, where the heavy mediator was integrated out. The appropriate EFT for this purpose is known as non-relativistic QED (NRQED) [68]. This EFT is often used to study nuclear effects in atoms, but here we will use it to describe our heavy mediator. In this theory the effects of the heavy mediator are captured by four-fermion interactions suppressed by powers of m φ . The reader who is familiar with NRQED may be surprised by the m −3 φ because the theory contains operators that are suppressed by m −2 φ and by m −4 φ but none that are suppressed by a cubic power [69]. This however can be resolved within NRQED as we now show, since the m −3 φ is simply a 1-loop correction to the m −2 φ Wilson coefficient 5 . To demonstrate the m −3 φ scaling one needs to do a 1-loop matching between the full Yukawa theory and the EFT. This calculation amounts to calculating the LO QED correction to the tree level calculation in the Yukawa theory.
To perform the matching we calculate a physical quantity twice, once in the full theory and once in the EFT, and then set the Wilson coefficient in the EFT to get a similar result order by order in perturbation theory. For the full theory we take non-relativistic QED with an additional Yukawa interaction. For simplicity we take the nucleus to be infinitely heavy and of zero size. The nuclear charge will again be set to Z eff , in accordance with the approximation taken above for electron screening effects. The Yukawa propagator in momentum space is (q 2 + m 2 φ ) −1 and the rest of the Feynman rules are shown in [70].
The EFT we will match onto is NRQED, again with an infinitely heavy nucleus, but with the inclusion of one contact interaction between the nucleus, N , and the electron C(ψ * e ψ e )(ψ * N ψ N )/m 2 φ , where C is a Wilson coefficient. Since we are only interested in the scaling behavior we will not dwell on precise numerical coefficients, but will focus on the parametric scaling. A quantity that is simple to compute for this matching calculation is the two-electron correlation function setting external momenta to zero, thus avoiding complications of bound states.
We begin by matching at tree level. In the full theory, to LO in the NP interaction, the tree-level contribution to the 2-point function is shown in Fig. S3 and is trivially α NP /m 2 φ . In the EFT, a similarly simple calculation produces C/m 2 φ for the two-point function. As expected, at tree level the matching gives C tree = α NP (S13) Still at LO in NP, we can calculate to an additional order in the QED coupling α which arrises at one loop and is also shown in Fig. S3. Since Coulomb scattering diverges as the incoming velocity goes to zero, we will not be surprised to encounter an IR divergence in this calculation. However, as expected, the IR divergence has an identical structure in the full and effective theories and will thus not affect the value of the Wilson coefficient at order α. In the full theory we find the two point function Full theory: where the three terms in the integrand are the Yukawa, electron and Coulomb propagator, respectively, and we used α = e 2 /4π and E=0. The coulomb potential is IR-regulated by λ. The result in Eq. (S14) has been expanded at large m φ . The second term in Eq. (S14) is already reminiscent of the m −3 φ term in Eq. (S12). Repeating the O(α) computation in the EFT is straightforward, Z eff e 2 q 2 + λ 2 ∼ − 2CZ eff αm e λm 2 φ . (S15) As expected, using the tree-level value for C, Eq. (S13), the IR divergences in the full and effective theories match at order α. However, the second term in Equation (S14), leads to a finite correction to the Wilson coefficient at O(α), We thus find that the 1-loop correction to the m −2 φ Wilson coefficient shown in Fig. S3 leads to an effective operator which scales as Z eff /(a 0 m φ ) 3 , in qualitative agreement with our position-space calculation. The two calculations are related. The position space calculation makes use of the hydrogen wavefunction and thus in some sense re-sums diagrams with an arbitrary number of photon exchanges between the nucleus and the electron before and after the scalar line. The EFT calculation isolates the first of these corrections. In both calculations we have accounted for multi-electron effects by taking the nuclear charge to be an effective one, with the potential being Coulombic otherwise. Again, we find that the mismatch between the Z eff 's does not allow this term to be absorbed in the FS. We are, however, not surprised that the numerical coefficients of the m −3 φ terms in Eqs. (S12) and (S16) do not agree since the former captures bound state dynamics, properly resumming IR effects that are dominant at the Bohr radius.