Probing new light force-mediators by isotope shift spectroscopy

TL;DR

This paper proposes using isotope shift spectroscopy to detect new light force-mediators, developing a formalism to interpret data as bounds on new physics, and demonstrating potential for high sensitivity in future experiments.

Contribution

It introduces a novel formalism for interpreting isotope shift data as bounds on new light bosons, independent of nuclear effects, and applies it to existing and future experimental data.

Findings

Current Ca+ data sets bounds on new bosons.

Projected future measurements could improve sensitivity by five orders of magnitude.

Potential to detect bosons in the 10 keV to 10 MeV mass range.

Abstract

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.

Tables

Table 1. Table 1: The 95 % CL bounds on y e ​ y n subscript 𝑦 𝑒 subscript 𝑦 𝑛 y_{e}y_{n} for a massless mediator ϕ italic-ϕ \phi from Ca + data  PhysRevLett.115.053003 and 95 % CL projections for Ca + , Sr + , Sr/Sr + and Yb + assuming on error of σ i = 1 subscript 𝜎 𝑖 1 \sigma_{i}=1\, Hz.

	transition 1	transition 2	accuracy	$y_e{y}_n$ bound
	[nm]	[nm]	${\sigma }_i$	($m_{\varphi }=0$)
Ca+	397.0	866.5	$0.1\mathrm{MHz}$	$2\cdot {10}^{-9}$
Ca+	729.3	732.6	$1\mathrm{Hz}$	$2\cdot {10}^{-14}$
Sr+	674.0	687.0	$1\mathrm{Hz}$	$2\cdot {10}^{-13}$
Sr/Sr+	698.4	674.0	$1\mathrm{Hz}$	$3\cdot {10}^{-15}$
Yb+	435.5	466.9	$1\mathrm{Hz}$	$2\cdot {10}^{-15}$

Equations

$${\nu }_i^{A{A}^{\prime }}\equiv {\nu }_i^A-{\nu }_i^{A^{\prime }}.$$

$${\nu }_i^{A{A}^{\prime }}=K_i{\mu }_{A{A}^{\prime }}+F_i\delta ⟨r^2{⟩}_{A{A}^{\prime }}+\dots ,$$

$$m{\nu }_2^{A{A}^{\prime }}=K_{21}+F_{21}m{\nu }_1^{A{A}^{\prime }},$$

$${\overrightarrow{m\nu }}_i\equiv \left(m{\nu }_i^{A{A}_1^{\prime }},m{\nu }_i^{A{A}_2^{\prime }},m{\nu }_i^{A{A}_3^{\prime }}\right).$$

$${\overrightarrow{m\nu }}_i=K_i\overrightarrow{m\mu }+F_i\overrightarrow{m\delta ⟨r^2⟩}.$$

$$\mathrm{NL}=\frac{1}{2}\left|({\overrightarrow{m\nu }}_1\times {\overrightarrow{m\nu }}_2)\cdot \overrightarrow{m\mu }\right|.$$

$${\nu }_i^{A{A}^{\prime }}=K_i{\mu }_{A{A}^{\prime }}+F_i\delta ⟨r^2{⟩}_{A{A}^{\prime }}+{\alpha }_{\mathrm{NP}}{X}_i{\gamma }_{A{A}^{\prime }},$$

$${\overrightarrow{m\nu }}_2=K_{21}\overrightarrow{m\mu }+F_{21}{\overrightarrow{m\nu }}_1+{\alpha }_{\mathrm{NP}}\vec{h}{X}_1\left(X_{21}-F_{21}\right),$$

$${\alpha }_{\mathrm{NP}}$$

$${\mathrm{NL}}_{\mathrm{NP}}=\frac{{\alpha }_{\mathrm{NP}}}{2}(\overrightarrow{m\mu }\times \vec{h})\cdot (X_1{\overrightarrow{m\nu }}_2-X_2{\overrightarrow{m\nu }}_1).$$

$$X_i=\frac{1}{A-Z}{\frac{d{ϵ}_{ab}}{d{\alpha }_{\mathrm{NP}}}|}_{{\alpha }_{\mathrm{NP}}=0},$$

$$X_i=\int d^{3}r\frac{e^{-m_{\varphi }r}}{r}\left[|{\Psi }_b(r){|}^2-|{\Psi }_a(r){|}^2\right],$$

$${\overrightarrow{m\nu }}_i$$

$$\overrightarrow{m\delta ⟨r^2⟩}$$

$$\overrightarrow{m\mu }$$

$$\alpha Z\to \alpha Z+{\alpha }_{\mathrm{NP}}(A-Z).$$

$$E_a^A=(\alpha Z_{\mathrm{eff}}^a{)}^2{I}_a^A,$$

$$\alpha Z_{\mathrm{eff}}^a\to \alpha Z_{\mathrm{eff}}^a+{\alpha }_{\mathrm{NP}}(A-Z-{\sigma }_a^{\prime }),$$

$${\nu }_i^{A{A}^{\prime }}=$$

$$=$$

$$\displaystyle\nu^{AA^{\prime}}_{i}\approx\nu^{AA^{\prime}}_{i}\big{|}_{\alpha_{\rm NP}=0}+\frac{2\alpha_{\rm NP}}{\alpha}\left[A\left(\frac{E^{A}_{b}}{Z_{\rm eff}^{b}}-\frac{E^{A}_{a}}{Z^{a}_{\rm eff}}\right)-(A\to A^{\prime})\right]\,.$$

$${X_i|}_{m_{\varphi }=0}\approx 2{\alpha }^{-1}\left(\frac{E_b}{Z_{\mathrm{eff}}^b}-\frac{E_a}{Z_{\mathrm{eff}}^a}\right).$$

$$[{\sigma }_{{\alpha }_{\mathrm{NP}}}{]}_{\text{proj}}\sim \frac{\sqrt{{\sigma }_2^2+{\sigma }_1^2{F}_{21}^2}}{\left(X_2-X_1{F}_{21}\right)}\frac{A}{\Delta A_j^{\min }\Delta A_j^{\max }},$$

$$\Delta E_a={\alpha }_{\mathrm{NP}}\int d^{3}r\frac{e^{-m_{\varphi }r}}{r}{\left|{\Psi }_a(r)\right|}^2.$$

$$\Delta E_{n\mathrm{S}}={\alpha }_{\mathrm{NP}}\frac{{\left|{\psi }_{n\mathrm{S}}(0)\right|}^2}{m_{\varphi }^2}\left(1-\frac{4Z_{\mathrm{eff}}^{(n\mathrm{S})}\alpha m_e}{m_{\varphi }}+\dots \right),$$

$$C_{\mathrm{tree}}={\alpha }_{\mathrm{NP}}$$

$$\text{mbox}Fulltheory:D_{ee}^{(1)}=\int \frac{d^{3}q}{(2\pi {)}^3}\frac{{\alpha }_{\mathrm{NP}}}{q^2+m_{\varphi }^2}\frac{1}{E-\frac{q^2}{2m_e}}\frac{Z_{\mathrm{eff}}{e}^2}{q^2+{\lambda }^2}\sim -\frac{2{\alpha }_{\mathrm{NP}}{Z}_{\mathrm{eff}}\alpha m_e}{\lambda m_{\varphi }^2}+\frac{2{\alpha }_{\mathrm{NP}}{Z}_{\mathrm{eff}}\alpha m_e}{m_{\varphi }^3}+\dots$$

$$\text{mbox}EFT:D_{ee}^{(1)}=\int \frac{d^{3}q}{(2\pi {)}^3}\frac{C}{m_{\varphi }^2}\frac{1}{E-\frac{q^2}{2m_e}}\frac{Z_{\mathrm{eff}}{e}^2}{q^2+{\lambda }^2}\sim -\frac{2C{Z}_{\mathrm{eff}}\alpha m_e}{\lambda m_{\varphi }^2}.$$

$$\Delta C=\frac{2{\alpha }_{\mathrm{NP}}{Z}_{\mathrm{eff}}\alpha m_e}{m_{\varphi }}=\frac{2{\alpha }_{\mathrm{NP}}{Z}_{\mathrm{eff}}}{a_0{m}_{\varphi }}.$$

Full text

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.

††preprint: DESY 17-055, FERMILAB-PUB-17-077-T, LAPTh-009/17, MIT-CTP-4898

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 matter-antimatter 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.: Baron:2013eja ; PhysRevLett.100.120801 ; bouchendira2011new ; wood1997measurement ; guena2005measurement , for a review see Budker:2017ReviewInPrep ; Karshenboim20051 ; Ginges:2003qt ) 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 Kim:1986ax ; Graham:2015cka ; Gupta:2015uea ; Flacke:2016szy ; Feng:1997tn ; Wilczek:1982rv ; Gelmini:1982zz . 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 King:63 of isotope shifts of two narrow transitions Delaunay:2016brc . 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-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 $\nu$ comparing the isotopes $A$ and $A^{\prime}$ is the IS,

[TABLE]

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

[TABLE]

where the two terms represent MS and FS respectively King:63 ; King:13 . We define $\mu_{AA^{\prime}}\equiv m^{-1}_{A}-m^{-1}_{A^{\prime}}$, where $m_{A}$ and $m_{A^{\prime}}$ are the masses of isotopes $A$ and $A^{\prime}$.

The quantity $\delta\langle r^{2}\rangle_{AA^{\prime}}$ is dominated by the difference in the mean squared charge radii of the two nuclei but can include other contact interactions. Both $\mu_{AA^{\prime}}$ and $\delta\langle r^{2}\rangle_{AA^{\prime}}$ are purely nuclear quantities that do not depend on the electronic transition $i$. Note, however, that $\mu_{AA^{\prime}}$ is known with high precision, whereas $\delta\langle r^{2}\rangle_{AA^{\prime}}$ 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 $\delta\langle r^{2}\rangle_{AA^{\prime}}$ giving a relation between the isotope shifts $\nu_{1}^{AA^{\prime}}$ and $\nu_{2}^{AA^{\prime}}$. In terms of the modified IS111Below we will adopt the notation of adding an $m$ to “modified” (i.e. normalized by $\mu_{AA^{\prime}}$) quantities, such as $m\delta\langle r^{2}\rangle_{AA^{\prime}}\equiv\delta\langle r^{2}\rangle_{AA^{\prime}}/\mu_{AA^{\prime}}$., $m\nu_{i}^{AA^{\prime}}\equiv\nu_{i}^{AA^{\prime}}/\mu_{AA^{\prime}}$, this relation is,

[TABLE]

with $F_{21}\equiv F_{2}/F_{1}$, and $K_{21}\equiv K_{2}-F_{21}K_{1}$.

Equation (3) reveals a linear relation between $m\nu_{1}$ and $m\nu_{2}$, giving rise to a straight line in the so-called King plot of $m\nu_{2}$ vs $m\nu_{1}$  King:63 . 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

[TABLE]

The nuclear parameters of field and mass shift, $\mu_{AA^{\prime}}$ and $\delta\langle r^{2}\rangle_{AA^{\prime}}$ can also be written as vectors $\overrightarrow{m\mu}$ and $\overrightarrow{m\delta\langle r^{2}\rangle}$ in the same space (notice that $\overrightarrow{m\mu}\equiv\left({1,1,1}\right)$) and hence Eq. (2) becomes

[TABLE]

In this language LO factorization implies the following qualitative statement: any vector of reduced isotope shifts, $\overrightarrow{m\nu}_{i}$, must lie in the plane that is defined by $\overrightarrow{m\mu}$ and $\overrightarrow{m\delta\langle r^{2}\rangle}$, as illustrated in the cartoon in the left panel of Fig. S1.

Note that, because the direction of $\overrightarrow{m\delta\langle r^{2}\rangle}$ in this space is uncertain, theory does not tell us in which direction this plane is oriented. However, by measuring two IS vectors, $\overrightarrow{m\nu}_{1}$ and $\overrightarrow{m\nu}_{2}$, we can test this statement by asking whether the three vectors $\overrightarrow{m\nu}_{1}$, $\overrightarrow{m\nu}_{2}$, and $\overrightarrow{m\mu}$ are coplanar. The coplanarity of these vectors corresponds to King linearity as we can see by rewriting Eq. (3) in vectorial form $\overrightarrow{m\nu}_{2}=K_{21}\overrightarrow{m\mu}\!+\!F_{21}\overrightarrow{m\nu}_{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 $\overrightarrow{m\nu}_{1}$ and $\overrightarrow{m\nu}_{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,

[TABLE]

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 $\overrightarrow{m\nu}_{1,2}$ and $\overrightarrow{m\mu}$. A given data set is considered linear if NL is smaller than its first-order propagated error $\sigma_{\rm NL}=\sqrt{\Sigma_{k}(\partial\textrm{NL}/\partial{\cal O}_{k})^{2}\sigma_{k}^{2}}$ where the sum runs over all measured observables ${\cal O}_{k}$ (modified frequency shifts and isotope masses) with standard deviations $\sigma_{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),

[TABLE]

namely $X_{i}$ depends on the form of the new potential and on the electronic transition, while $\gamma_{AA^{\prime}}$ depends only on the nuclear properties. The parameter $\alpha_{\rm 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}\propto F_{i}$. In this case the NP term can be absorbed by redefining $\delta\langle r^{2}\rangle_{AA^{\prime}}$. Also, if the new physics couples to electrons and nuclei according to their electric charge (such as the case of dark-photon Holdom:1985ag ), $\gamma_{AA^{\prime}}=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

[TABLE]

where $\vec{h}$ is the NP vector in reduced frequency units, that is $h_{AA^{\prime}}\equiv\gamma_{AA^{\prime}}/\mu_{AA^{\prime}}$ and $X_{21}\equiv 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}\!\neq\!F_{21}$; (ii) $\vec{h}$ is not aligned with any linear combination of $\overrightarrow{m\mu}$, $\overrightarrow{m\nu}_{1}$ or $\overrightarrow{m\nu}_{2}$.

By solving the set of equations (7) one finds an expression for $\alpha_{\rm NP}$ that is needed to yield a particular dataset $\left\{\overrightarrow{m\nu}_{1},\overrightarrow{m\nu}_{2},\overrightarrow{m\mu}\right\}$,

[TABLE]

assuming NP is the dominant contribution to nonlinearity. If linearity holds then $\alpha_{\rm NP}\lesssim\sigma_{\alpha_{\rm NP}}=\sqrt{\Sigma_{k}(\partial\alpha_{\rm NP}/\partial{\cal O}_{k})^{2}\sigma_{k}^{2}}$. Hence, the sensitivity to probe $\alpha_{\rm NP}$ is lost in the limit where the denominator in Eq. (9) vanishes, because the NP contribution to nonlinearity is

[TABLE]

It is straightforward to check that this happens under the conditions specified below Eq. (8).

The presented method of limiting $\alpha_{\rm NP}$, Eq. (9), contains theory input only in $X_{i}$ and $h_{AA^{\prime}}$ which describe how NP affects the IS. The SM contribution in the factorized limit is fully parametrized by the observables $\vec{\nu}_{i}$ and $\vec{\mu}$. The form of $h_{AA^{\prime}}$ 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^{\prime}}\propto AA^{\prime}$ Haber:1978jt ; Delaunay:2016brc . The atomic transition-dependent factors $X_{1,2}$ can be reasonably calculated by a many-body 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 SmExp:1981 ). Nevertheless, some level of nonlinearity is expected to arise from SM higher-order contributions 0022-3700-15-7-009 ; PhysRev.188.1916 ; 0022-3700-20-15-015 ; PhysRevA.31.2038 . 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 ${\cal O}(10^{-3}-10^{-4})$ in light atoms Puchalski2010 . Likewise, such contributions in heavier elements with $Z=20-87$ 0022-3700-20-15-015 , but only for $S\to 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 optical-clock transitions in laser-cooled atoms or ions JunYe ; NIST ; NRC ; Innsbruck ; PTB ; PhysRevLett.116.063001 . Second, since the hyperfine interaction of electrons with the nucleus is a source for King nonlinearity 0022-3700-15-7-009 , 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_{\phi}$ 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_{\phi}(r)=-\alpha_{\rm NP}(A-Z)e^{-m_{\phi}r}/r$, where $\alpha_{\rm NP}=(-1)^{s}y_{e}y_{n}/4\pi$. 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) dzuba96pra . For the single-valence electron ions Ca+ and Sr+, we create an operator $\hat{\Sigma}$ (see for example Berengut:2003zg ) 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\to b$) can be expressed

[TABLE]

where $\epsilon_{ab}$ is the difference of the energy levels of the states $a,~{}b$, evaluated as a function of $\alpha_{\rm NP}$ and the derivative is taken numerically at $\alpha_{\rm NP}=0$.

For neutral Sr, which has two valence electrons above closed shells, we use the CI+MBPT method as described in Berengut:2005bz . 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^{N_{e}-2}$ potential where $N_{e}$ is the total number of electrons). In this potential we generate a set of B splines PhysRevA.37.307 ; PhysRevLett.93.130405 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\rightarrow 4f^{13}\,6s^{2}$. For this ion we use the particle-hole CI+MBPT method in the $V^{N_{e}-1}$ potential berengut16pra which has previously been used for Hg+.

The many-body calculations can be cross-checked by perturbation theory, which yields

[TABLE]

where $|\Psi(r)|^{2}$ is the electron-density evaluated in the absence of NP, and $h_{AA^{\prime}}=AA^{\prime}\,$amu for the NP contribution in Eq. (7). As a cross-check of our many-body calculation, we use GRASP2K Jonsson20132197 to evaluate $|\Psi(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, $r_{N}\sim\,A^{1/3}\times(200\,$MeV$)^{-1}$. Here $a_{0}\approx(4\,$keV$)^{-1}$ is the Bohr radius and $n_{e}$ is the ionization number. For $m_{\phi}\lesssim(1+n_{e})/a_{0}$, the “massless limit”, the interaction range is larger than the atomic size and $V_{\phi}\propto 1/r\,$ so that $X_{i}$ becomes independent of $m_{\phi}$. For intermediate masses, $(1+n_{e})/a_{0}\lesssim m_{\phi}\lesssim 1/r_{N}$, the interaction range is within the size of the electron wavefunction, and the potential $V_{\phi}\propto e^{-m_{\phi}r}/r$ is mass-dependent. Hence, detailed knowledge of the electronic wavefunctions is necessary to evaluate the effect of NP. In the heavy mass limit, $m_{\phi}\gtrsim 1/r_{N}$, the interaction range is shorter than the nuclear radius and $V_{\phi}\propto\delta(r)/(m_{\phi}^{2}r^{2})$. In this limit, the NP and nuclear charge-radius effects are approximately aligned since $X_{i}\propto F_{i}\propto|\Psi_{b}(0)|^{2}-|\Psi_{a}(0)|^{2}$. This results in a suppressed sensitivity for new physics which scales as $(X_{21}-F_{21})\to 0$, see Eq. (3), and Delaunay:2016brc .

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 $\alpha$, see Appendix II. We do not estimate the bounds on $\alpha_{\rm 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 ${\cal 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. 1.

V.1 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_{\phi}\lesssim 10\,$keV, the bound is essentially independent of $m_{\phi}$. 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_{\phi}$ dependence of the $y_{n}y_{e}$ bound at that limit. However, the ratio $F_{21}^{\rm th}/X_{21}$ ($F^{\rm 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_{\phi}$ increases, we rescale the $y_{e}y_{n}$ bound by $(1-F^{\rm exp}_{21}/X_{21})/(1-F_{21}^{\rm th}/X_{21})$, where $F_{21}^{\rm exp}$ 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_{\phi}>Z\alpha m_{e}$ the limits get weak, and the sensitivity decreases approximately as $m^{-3}_{\phi}$ for large masses. In Appendix IV we give two heuristic arguments that obtain this asymptotic scaling of our loss of sensitivity222We thank Richard Hill and Clara Peset for an enlightening conversation: the first is based on approximating Eq. (12); and the second is based on a non-relativistic QED (NRQED) effective theory approach.

For current bounds, we consider Ca+ ($Z=20$). There are five zero-nuclear-spin, stable or long-lived isotopes with $A=40,42,44,46,48$. Refs. PhysRevLett.115.053003 ; Shi2016 reported IS measurements for three isotope pairs ($A=42,44,48$ relative to 40) in three dipole-allowed transitions in Ca+ at wavelengths of 397.0 nm ($S\to P$), 866.5 nm ($D\to P$) and 854 nm ($D\to P$, not used here) with an uncertainty of ${\cal O}(100)\,\textrm{kHz}$.

Among the non-IS experiments that probe $y_{e}y_{n}-m_{\phi}$ 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$\times 100\,$eV, fifth force experiments Bordag:2001qi ; bordag2009advances 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 Harnik:2012ni .

Furthermore, separately $y_{n}$ is constrained by various neutron scattering experiments Leeb:1992qf ; Nesvizhevsky:2007by ; Barbieri:1975xy and $y_{e}$ by the anomalous magnetic moment of the electron $(g-2)_{e}$ Olive:2016xmw ; Hanneke:2010au and by electron beam-dump experiments for $m_{\phi}>1$ MeV.

Both $y_{e}$ and $y_{n}$ are also severely constrained by globular cluster energy loss for masses $m_{\phi}\lesssim 10$ keV Raffelt:2012sp ; Yao:2006px ; Grifols:1986fc ; Grifols:1988fv down to $y_{e}y_{n}<10^{-25}$ and $y_{e}$ by sun cooling Redondo:2008aa ; Jaeckel:2010ni . Couplings to nucleons in the $10^{-10}-10^{-7}$ range for $m_{\phi}\lesssim 100$ MeV may be also excluded by energy loss in the core of SN1987A Raffelt:2012sp ; Blum:2016afe . However, such astrophysical bounds might be avoided in certain models such as chameleon, see Burrage:2016bwy 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.

V.2 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 fine-structure 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_{\phi}$ and summarized in the lower part of Tab. 1. 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_{\phi}\geq 10\,\textrm{keV}$ region in parameter space. For lower $m_{\phi}$ 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 Burrage:2016bwy . Thus, an independent laboratory bound in this low-mass region is nevertheless worthwhile. For $m_{\phi}\sim$ a few MeV the projections of Ca+ ($S\to D$ transitions) and Sr+ are comparable to $y_{e}y_{n}$ from neutron scattering Leeb:1992qf and $(g-2)_{e}$. Since neutron experiments are affected by uncertainties Wissmann:1998ta ; Nesvizhevsky:2007by ; Antoniadis:2011zza ; TuckerSmith:2010ra 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\,\textrm{keV}$ Leeb:1992qf 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 $\sim 10\,\textrm{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 Feng:2016jff ; Feng:2016ysn can be probed by future IS measurements of Yb+ at the 1 Hz level.

Acknowledgment

We thank J. Feng, O. Firstenberg, I. Galon, R. Hill, S. G. Karshenboim, K. Pachucki, C. Peset, M. Pospelov, P. O. Schmidt and T. Tait for fruitful discussions. CD and YS are particularly thankful to J. Ekman and P. J$\ddot{{\rm o}}$nsson for their help regarding the GRASP2K program. The work of DB is supported by the DFG Reinhart Koselleck project, and the ERC (Dark-OST Advanced Project). The work of CD is suported by the “Investissements d’avenir, Labex ENIGMASS”. The work of CG is supported by the European Commission through the Marie Curie Career Integration Grant 631962 and by the Helmholtz Association through the recruitment initiative programme. The work of RO is supported by grants from the ISF, I-Core, ERC, Minerva and the Crown photonics center. The work of GP is supported by grants from the BSF, ERC, ISF, Minerva, and the Weizmann-UK Making Connections Programme. The work of YS is supported by the U.S. Department of Energy (DOE) under grant contract numbers DE-SC-00012567 and DE-SC-00015476.

I Visualizing the vector space

In the main text we define the following vectors in the $A^{\prime}$ vector space

[TABLE]

As long as $\overrightarrow{m\nu}_{1,2}$ are spanned by $\overrightarrow{m\mu}$ and $\overrightarrow{m\delta\langle r^{2}\rangle}$, the resulting King plot will be linear. In Fig. S1, we illustrate the vector space of the various components related to isotope shifts that leads to the nonlinearites. The NP contribution to IS, $\alpha_{\rm NP}X_{i}\vec{h}$, may lift the IS vectors from the $(\overrightarrow{m\mu},\overrightarrow{m\delta\langle r^{2}\rangle})$ plane, resulting in a nonlinear King plot. Fig. S2 illustrates a nonlinear King plot, where the area of the triangle corresponds to the NL of Eq. (6).

II Derivation of $X_{i}$ in the $m_{\phi}\to 0$ limit

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_{\phi}\ll(1+n_{e})/a_{0}\sim\mathcal{O}($few keV). Since in this limit the effective potential is Coulomb-like, $V_{\phi}(r)\simeq(A-Z)\alpha_{\rm NP}/r$, $X_{i}$ can be simply estimated through a shift of the fine-structure constant $\alpha$, without a detailed calculation of the electronic wavefunctions. At fundamental level, the Coulomb potential is modified by the shift

[TABLE]

In the absence of NP, the binding energy of the atomic level $a$ for isotope $A$ scales as

[TABLE]

where $Z^{a}_{\rm eff}\equiv Z-\sigma_{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 ${\cal O}(\alpha Z_{\rm eff}^{a})^{4}$ corrections from the fine-structure). The constant $\sigma_{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

[TABLE]

where the constant $\sigma^{\prime}_{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\to b$ is (in natural units, i.e. the reduced Planck constant is set to $\hbar=1$)

[TABLE]

and expanding to leading order in $\alpha_{\rm NP}\ll\alpha$ yields

[TABLE]

Since IS are typically orders of magnitude smaller than the transition frequencies (with the exceptions of very degenerate states such as in dysprosium PhysRevA.50.132 ; Leefer:2014tga ), we can take $E^{A^{\prime}}\approx E^{A}$ in the NP contribution above and matching to Eq. (7) we find $\gamma_{AA^{\prime}}=A-A^{\prime}$ and

[TABLE]

By a combination with Eq. (9), the above expression provides a reasonable estimate of the constraint on $\alpha_{\rm 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 ${\cal 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 $\sigma_{i}$ of $\nu_{i}^{AA^{\prime}}$, the only missing information is how much the new physics vector, $\vec{h}$, points towards the linearity plane. While a precise determination of the vector component requires data, the projection of the NP vector $h_{AA^{\prime}}=AA^{\prime}\,$amu along the mass shift direction is easily obtained without the need of experimental input. Therefore, a best-case projection333The actual bound obtained by data will be always weaker, $[\sigma\alpha_{\rm NP}]_{\text{proj}}\leq[\sigma\alpha_{\rm NP}]_{\text{data}}$ since the projection neglects the alignment with the FS that will weaken the bound. $[\sigma_{\alpha_{\rm NP}}]_{\text{proj}}$ can be obtained by neglecting the possible additional alignment of NP with nuclear effects,

[TABLE]

where $\Delta A_{j}^{\rm min(max)}\equiv{\rm min\,(max)}[A-A_{j}]$ and $\sigma_{i}$ is the assumed standard deviation of IS measurements in transition $i$. Note that the sensitivity to $\alpha_{\rm NP}$ is weaker by a factor of $A/\Delta A_{j}^{\rm min}$ than the naive expectation of Refs. Delaunay:2016brc ; Frugiuele:2016rii ; Delaunay:2016zmu . 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/(\Delta A_{j}^{\rm min}\Delta A_{j}^{\rm max})$ 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\alpha m_{e})^{-1}=a_{0}/Z$. We thus expect our method to start losing power for mediators heavier than $Z\alpha m_{e}$ as can be seen in Fig. 1.

It is instructive to investigate like what power of $m_{\phi}/(Z\alpha 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}\propto 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_{\phi}\rightarrow\infty$. 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_{\phi}>Z\alpha m_{e}$ agrees with our theoretical expectation.

The asymptotic behavior of the limits in Fig. 1 exhibits an approximate $m_{\phi}^{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

[TABLE]

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)}_{\mathrm{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 $n$S state, with ($n=1,2,\ldots$), 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

[TABLE]

where we assumed that $m_{\phi}\gg Z\alpha m_{e}$ and expanded in powers of $1/m_{\phi}$, 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 $\delta^{3}(r)/m_{\phi}^{2}$. This term will thus lead to isotope shifts in the $a\to b$ transition that are proportional to $|\Psi_{a}(0)|^{2}-|\Psi_{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_{\mathrm{eff}}$ leading to an IS proportional to $Z_{\mathrm{eff}}^{(a)}|\Psi_{a}(0)|^{2}-Z_{\mathrm{eff}}^{(b)}|\Psi_{b}(0)|^{2}$. This term will not cancel in $(X_{2}-X_{1}F_{21})$, and since it scales as $m_{\phi}^{-3}$ we would expect the bounds to lose power as $m_{\phi}^{3}$ at high mass444It is interesting to note, that for a single electron atom, where $Z_{\mathrm{eff}}=Z$ for all states, the isotope shift again scales as $|\Psi_{a}(0)|^{2}-|\Psi_{b}(0)|^{2}$. In this case the bound will scale as $m_{\phi}^{4}$ from the higher order terms that were dropped in Eq. (S12)..

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) Caswell:1985ui . 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_{\phi}$. The reader who is familiar with NRQED may be surprised by the $m_{\phi}^{-3}$ because the theory contains operators that are suppressed by $m_{\phi}^{-2}$ and by $m_{\phi}^{-4}$ but none that are suppressed by a cubic power Hill:2012rh . This however can be resolved within NRQED as we now show, since the $m_{\phi}^{-3}$ is simply a 1-loop correction to the $m_{\phi}^{-2}$ Wilson coefficient555We thank Richard Hill for an enlightening conversation..

To demonstrate the $m_{\phi}^{-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_{\mathrm{eff}}$, in accordance with the approximation taken above for electron screening effects. The Yukawa propagator in momentum space is $(q^{2}+m_{\phi}^{2})^{-1}$ and the rest of the Feynman rules are shown in Kinoshita:1995mt .

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(\psi_{e}^{*}\psi_{e})(\psi_{N}^{*}\psi_{N})/m_{\phi}^{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 $\alpha_{\rm NP}/m_{\phi}^{2}$.

In the EFT, a similarly simple calculation produces $C/m_{\phi}^{2}$ for the two-point function. As expected, at tree level the matching gives

[TABLE]

Still at LO in NP, we can calculate to an additional order in the QED coupling $\alpha$ 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 $\alpha$. In the full theory we find the two point function

[TABLE]

where the three terms in the integrand are the Yukawa, electron and Coulomb propagator, respectively, and we used $\alpha=e^{2}/4\pi$ and E=0. The coulomb potential is IR-regulated by $\lambda$. The result in Eq. (S14) has been expanded at large $m_{\phi}$. The second term in Eq. (S14) is already reminiscent of the $m_{\phi}^{-3}$ term in Eq. (S12). Repeating the $O(\alpha)$ computation in the EFT is straightforward,

[TABLE]

As expected, using the tree-level value for $C$, Eq. (S13), the IR divergences in the full and effective theories match at order $\alpha$. However, the second term in Equation (S14), leads to a finite correction to the Wilson coefficient at $O(\alpha)$,

[TABLE]

We thus find that the 1-loop correction to the $m_{\phi}^{-2}$ Wilson coefficient shown in Fig. S3 leads to an effective operator which scales as $Z_{\mathrm{eff}}/(a_{0}{m_{\phi}})^{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_{\mathrm{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_{\phi}^{-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.