Probing new physics using Rydberg states of atomic hydrogen

We consider the role of high-lying Rydberg states of simple atomic systems such as H in setting constraints on physics beyond the Standard Model. We obtain highly accurate bound states energies for a hydrogen atom in the presence of an additional force carrier (the energy levels of the Hellmann potential). These results show that varying the size and shape of the Rydberg state by varying the quantum numbers provides a way to probe the range of new forces. By combining these results with the current state-of-the-art QED corrections, we determine a robust global constraint on new physics that includes all current spectroscopic data in hydrogen. Lastly we show that improved measurements that fully exploit modern cooling and trapping methods as well as higher-lying states could lead to a strong, statistically robust global constraint on new physics based on laboratory measurements only.


I. INTRODUCTION
Detailed measurements of atomic spectra were key to the discovery of quantum mechanics and the development of relativistic quantum electrodynamics (QED). Today, precision atomic spectroscopy underpins the SI system of units, provides the values of some fundamental constants, and enables precise tests of Standard Model calculations.
Looking for deviations between precise spectroscopic measurements and their Standard Model predictions thus provides a powerful way to set constraints on new physics [1]. One powerful approach looks for small effects that break symmetries such as parity (P violation) or timereversal (T-violation). Alternatively, one can compare experimental and theoretical transition frequencies. If additional force mediators (bosons) were present that coupled strongly enough to the nucleus and electrons, they would modify the frequency of spectral lines. Thus, by comparing experimentally measured spectra with theory the existence of new so-called fifth forces can be tested down to very small interaction strengths. In recent years extensions of the Standard Model, e.g., modified gravity [2,3], axions [4,5], new gauge boson [6,7], have been constrained in this way. In particular if the force mediator X is light, i.e., below 1 MeV in mass, and couples to partons and electrons, the limits obtained from atomic spectroscopy are many orders of magnitude stronger than from any other laboratory-based experiment, including high-energy collider experiments [6,8,9].
While modifications of the Standard Model through light bosons are predicted by various models, they arguably receive strong constraints from astrophysical * m.p.a.jones@durham.ac.uk † r.m.potvliege@durham.ac.uk ‡ michael.spannowsky@durham.ac.uk sources [10][11][12], e.g. the energy loss from the sun, globular clusters or supernovae. However, the need for independent laboratory-based experiments has been pointed out frequently -see, e.g., [13][14][15][16]. As an example, a prominent class of light-scalar models potentially related to modified gravity and dark energy, are chameleons [16][17][18]. Chameleons have a mass that depends on the energy density of their environment and thus can avoid being produced in stars, thereby avoiding astrophysical bounds.
One of the main uncertainties in the prediction of spectral lines arises due to the difficulty of solving the Schrödinger or Dirac equations for many interacting electrons. Even state-of-the-art calculations for species commonly used in atomic clocks only attain a fractional uncertainty of ∼ 10 −5 [19], which is 14 orders of magnitude lower than the current experimental precision. To circumvent this limitation, it has been proposed to look for new physics using the difference in spectral line positions between isotopes (isotope shifts) [5,20], rather than by direct comparison with theory. Although promising [21], the method is limited by the requirement that at least three stable isotopes with two suitable transitions exist for each element.
An alternative approach is to use light atomic species such as H or He for which full standard model predictions of line positions including QED corrections (Lamb shift) and weak interactions (Z boson exchange) are possible. Even here however, the complex structure of the nucleus, in particular the details of its charge distribution, limits the achievable accuracy. Spectroscopic data that does not strongly depend on the details of how the nucleus is modelled can thus help to improve the sensitivity on the presence of new forces.
In this paper, we explore how the precision spectroscopy of states with a high principal quantum number n (Rydberg states) might be used to set constraints on physics beyond the Standard Model. In principle such states offer several advantages that could be exploited in arXiv:1909.09194v1 [hep-ph] 19 Sep 2019 a search for new physics. Firstly, the overlap of Rydberg wave functions with the nucleus scales as n −3 , vastly reducing their sensitivity to nuclear effects. The radiative lifetime also scales as n −3 , meaning that narrow transitions from low-lying atomic states are available that span the UV to NIR wavelength range that amenable to precision laser spectroscopy. The n −2 scaling of the energy levels means that for each atomic species a large number of such transitions are available within a narrow spectral range. Lastly, there is a natural length scale associated with the atomic wave function that scales as n 2 . As we will show, being able to vary this length scale enables tests which are sensitive to the corresponding length scale associated with any new forces [8].
Here we take hydrogen as a model system in which to explore the use of Rydberg states in searches for new physics. Measurements with a fractional uncertainty of 10 ppt or better are already available for n up to 12. We calculate the (non-relativistic) spectrum of the combination of a hydrogenic Coulomb potential and a Yukawa potential arising from new physics to high accuracy. By combining the resulting energies with previously derived relativistic, QED and hyperfine corrections, we obtain predicted atomic transition frequencies that can be compared directly to experimental data to set a constraint on the strength of a new physics interaction. We consider in detail how uncertainties due to the Rydberg constant and the proton charge radius can be reduced or eliminated altogether, and show how a global statistical analysis can be used to derive robust atomic physics constraints. Lastly we develop proposals for future improved tests using atomic Rydberg spectroscopy in atomic hydrogen and other species.
The structure of the paper is as follows. In Section II we introduce the simplified model used to parameterize the effect of new physics. Calculations of the effect of new bosons on atomic energy levels are presented in section III. We assess the current experimental reach for new physics (NP) in Section IV. In Section V we discuss the impact of potential experimental and theoretical improvements on the uncertainty budget and in how far this can result in tighter constraints of new physics. We offer a summary and conclusions in Section VI.

II. PARAMETRISATION OF NEW PHYSICS
With the discovery of the Higgs boson [22,23], for the first time a seemingly elementary scalar sector was established in nature. Such a particle would mediate a new short-ranged force, the so-called Higgs boson force [24]. While the Higgs boson force is very difficult to measure in atom spectroscopy [20], many extensions of the Standard Model predict elementary scalar or vector particles with a very light mass. Examples include axions [4,5,25], modified-gravity models [2,26], millicharged particles [27][28][29], Higgs-portal models [30,31] and light Z [32,33]. To remain as model-independent as possible in parameterizing deformations from the Standard Model (SM), it has become standard practice to express new physics contributions in terms of so-called simplified models [34]. The idea is to add new degrees of freedom to the Standard Model Lagrangian without asking how such states arise from a UV complete theory. Thus, one can describe the dynamics and phenomenological implications of new degrees of freedom without making further assumptions on the UV theory from which they descend.
For example, if we assume a fifth force to be mediated through a novel spin-0 particle X 0 that couples to leptons and quarks with couplings g li and g qi respectively, we can augment the Standard Model Lagrangian L SM to Here i denotes the three flavor generations and l i and q i refer to the mass basis of the SM fermions. We note that the interactions of Eq. (1) could be straightforwardly extended to (axial)vector or pseudoscalar particles and to flavor off-diagonal interactions, e.g. g qijqi q j X 0 with i = j. Further, we should emphasize that the operators of Eq. (1) are gauge invariant only after electroweak symmetry breaking, which implies that the coefficients g fij must implicitly contain a factor v/Λ NP (v is the vacuum expectation value of the Higgs field and Λ NP is a new physics scale). However, this is only important for the interpretation of the observed limit we derive on g fi . Studying Rydberg states in hydrogen atoms, we will set a limit on the combined interaction g e g N , where g e and g N corresponds to the interaction of X 0 respectively with the electron and the nucleon. With the Lagrangian of Eq. (1), the interaction mediated by the NP boson X 0 between these two particles contributes an additional Yukawa potential V (r) to the Hamiltonian. Denoting by r the distance between the electron and the nucleon and by m X0 the mass of the particle, where s, an integer, is the spin of the force mediator (e.g., s = 0 for a scalar particle). Higher integer-spin mediators would also give rise to a Yukawa potential of this form. There is however a subtle difference in the sign of this potential between even and odd integer-spin force carriers. Lorentz invariance and the unitarity of the transition matrix element lead to an attractive (repulsive) force if g e g N > 0 (g e g N < 0) in the case of an even-spin mediator, and to an attractive (repulsive) force if g e g N < 0 (g e g N > 0) in the case of an odd-spin mediator. For example, as the charges for the Higgs boson (spin-0) and the graviton (spin-2) are the particles' masses, the Higgs force and gravity are both attractive. As we want to remain agnostic about the force carrier and the way it interacts with the nucleons and electrons, in the following we will allow both positive and negative values for g e g N . Finally, we note that an excellent recent review of this type of simplified model is provided in [1].

III. NEW PHYSICS LEVEL SHIFTS
The presence of the interaction potential V (r) would affect the atomic transition frequencies. Its effect can be evaluated perturbatively. To first order in V (r), and neglecting spin-orbit coupling and other relativistic corrections, the energy of a hydrogenic state of principal quantum number n, orbital angular momentum quantum number l and radial wave function R nl (r) is shifted by a quantity δE NP nl , with Since the interaction is spherically symmetric, the perturbation is diagonal in l and in the magnetic quantum number m, and δE NP nl does not depend on m. Taking into account V (r) to all orders, which we have done as a test of our numerical methods, confirms that second-and higher-order terms of the perturbation series are completely negligible for the couplings of interest, i.e. |g e g N | < 10 −11 .
The shift δE NP nl takes on a particularly simple form in the limit m X0 → 0: Since the virial theorem guarantees that where E n is the non-relativistic energy of the (n, l) states, α is the fine structure constant and Z is the number of protons in the nucleus. Moreover, (See Appendix A for the origin of the factor of 1/α and more generally for the conversion between natural units and atomic units.) Simple analytical forms can also be derived, e.g., for states with maximum orbital angular momentum (l = n − 1). However, in most cases δE NP nl is best evaluated numerically.
Various approaches to this problem have been considered over the years, as has the calculation of energy levels for a superposition of a Coulomb potential and a Yukawa potential (the Hellmann potential) [35][36][37][38][39][40][41][42][43]. The most accurate results reported to date are those of Ref. [40], in which the energies of the ground state and first few excited states were obtained to approximately 13 significant figures using a generalized pseudo-spectral method. Our approach to this problem is different and does not seem to have been used so far in this context: We expand the radial wave functions on a finite Laguerre basis of Sturmian functions S κ νl (r) [44], find the generalized eigenvectors of the matrix representing the unperturbed Hamiltonian in that basis, and use these to calculate the first order energy shift ∆E nl . Here with κ a positive parameter which can be chosen at will. These basis functions have already been used in this context, but in a different way [41]. Sturmian bases have proved to be convenient in precision calculations of properties of hydrogenic systems [45,46]. We obtain the eigenenergies and wave functions of the unperturbed Hamiltonian by solving the generalized eigenvalue problem where H 0 is the matrix representing the unperturbed nonrelativistic Hamiltonian of hydrogen in this basis and S is the overlap matrix of the basis functions (Sturmian functions are not mutually orthogonal). The corresponding matrix elements and the elements of the matrix V representing the Yukawa potential can be obtained in closed form using standard integrals and recursion formula [47]. Having the eigenvectors c, the energy shifts are then calculated as δE NP c = c T Vc. Since the functions {S κ νl (r), ν = 1, 2, . . .} form a complete set, the eigenvalues E and energy shifts δE NP c obtained with a basis of N of these functions (ν = 1, . . . , N ) converge variationally to the exact eigenenergies and exact energy shifts when N → ∞. We repeat the calculations for several different values of κ and different basis sizes so as to monitor the convergence of our results and the impact of numerical inaccuracies. With an appropriate choice of κ, and taking N up to 200, the calculated energy shifts converged to at least 8 significant figures [48]. Using the same method, but solving the generalized eigenvalue problem for the full Hamiltonian rather than the unperturbed Hamiltonian, we could also reproduce the results of Ref. [40] to the 14 significant figures given in that article.
The results of these calculations are summarized in Figs. 1, 2 and 3. These results, like all the other numerical results discussed in this paper, refer to the specific case of atomic hydrogen. We will therefore assume that Z = 1 from now on. Fig. 1 shows the general trends. The fractional shift is largest for light bosons, where the range of the Yukawa potential is comparable to or larger than the range of the atomic wave function. In agreement with Eq. (6), |δEnl NP /E n | |g e g N |/2πα for low masses. As m X0 increases, the shift decreases, but in a way that depends on the shape of the atomic wave function through both n and l. The effect of the Yukawa potential is largest at I. The range of the Yukawa potential (Λ), expressed as a multiple of the Bohr radius, and the principal quantum number nΛ for which this range is equal to that of the corresponding l = 0 state to the closest approximation possible, for three values of mX 0 , the mass of the NP particle.  To gain further insight, in Fig. 2, we investigate the relationship between the two characteristic length scales of the problem, i.e., the range of the Yukawa potential, Λ = 1/m X0 , and the range of the atomic wave function. The latter can be characterized by the expectation value nl|r|nl , which for l = 0 states is 3a 0 n 2 /2 where a 0 is the Bohr radius. We see that these two ranges are comparable for principal quantum numbers n ∼ n Λ , where n Λ is the integer closest to (2 Λ/3 a 0 ) 1/2 . Representative values of n Λ are given in Table I. The NP shift is accurately predicted by Eq. (6) for n n Λ and is much smaller than that limit for n n Λ . The fractional shift is plotted in Figs. 2(a) and (b), respectively against the ratio of these two characteristic lengths and against the boson mass, for a range of values of n and l. These curves show that for masses below ∼ 50 eV, the shift decreases with n but is essentially independent of l. Above this breakpoint, the shift decreases much more rapidly for dstates (l = 2) than for s-states (l = 0). This trend is even more marked for higher values of l (not shown in the figure). In fact, for states with l = n − 1 (which is the maximum value of the orbital angular momentum for the principal quantum number n), |δE NP nl /E n | decreases as fast as n −2n when n increases beyond n Λ .
In Figure 3, we fix the value of the fractional NP shift at |δE NP nl /E n | = 10 −12 , and show how the resulting constraint on the mass m X0 and the effective coupling g e g N depend on the quantum numbers n and l. Thus combining measurements for different values of n and l could provide additional information on the properties of the fifth-force carrier, i.e., its mass and its couplings to the electron and nuclei. FIG. 2. The NP shift, δE NP nl , divided by the non-relativistic energy of the state, En, for the states of atomic hydrogen with n = 5 (solid curves), n = 8 (dashed-dotted curves), n = 12 (dashed curves) or n = 26 (dotted curves) and l = 0 (black curves) or l = 2 (red curves), vs., (a) the range of the NP potential divided by the characteristic length scale of the atomic wave function, nl|r|nl , or (b) the mass of the NP particle. A value of gegN of 1 × 10 −12 is assumed.

IV. NP BOUNDS BASED ON CURRENT SPECTROSCOPIC DATA
In a nutshell, the existence of a new physics interaction could be brought to light by demonstrating a significant difference between the measured transition frequency for a transition from a state a to a state b, ∆ exp ba , and the corresponding prediction of the Standard Model, ∆ SM ba (or, better, by demonstrating such a difference for a set of transitions). Bounds on the strength of the new physics interaction can be set by finding the most positive and most negative values of g e g N for which ∆ exp ba is consistent with the theoretical value ∆ SM ba + ∆ NP ba with where h the Planck's constant. However, ∆ SM ba depends on the Rydberg constant R ∞ , whose value is primarily obtained by matching spectroscopic data to theory [49]. Setting bounds on g e g N makes it therefore necessary to evaluate ∆ SM ba with a value of R ∞ itself obtained with allowance made for the possibility of new physics shifts on the relevant atomic transitions. Frequency intervals have been both measured and calculated to a very high level of precision for transitions in hydrogen, deuterium and muonic hydrogen. However, a new physics interaction might couple an electron differently to a deuteron than to a proton, and couple a proton differently to a muon than to an electron. It is therefore prudent, when establishing such bounds, to use data pertaining to only one of these three systems rather than using mixed sets of data. We consider bounds based exclusively on hydrogen results in this paper.
∆ SM ba is the sum of a gross structure contribution ∆ g ba (as given by the elementary treatment based on the Schrödinger equation) and of various corrections arising from the Dirac equation, from QED effects and from the hyperfine coupling [49][50][51]. In terms of the Rydberg frequency, R = c R ∞ , where m r is the reduced mass of the atom and m e is the mass of the electron. It is convenient to factorize ∆ g ba into the product R∆ g ba , with The difference ∆ SM ba − ∆ g ba depends on R p , the charge radius of the proton, through a term roughly proportional to R 2 p [50,51]. We denote this term by R 2 p∆ ns ba , aggregate all the other corrections into a shift ∆ oc ba , and write The term ∆ oc ba includes fine structure and recoil corrections as well as QED and hyperfine shifts. Detailed work by a number of authors has yielded expressions for these corrections in terms of R, of R p and of a small number of fundamental constants determined from measurements in physical systems other than hydrogen. The values of R and R p recommended by the Committee on Data of the International Council for Science (CODATA) were co-determined by a global fit of the theory to a large set of data, including deuterium data [49]. Taking new physics shifts into account in a determination of R based entirely on hydrogen data thus involves a simultaneous redetermination of R p . Eq. (12) is a convenient starting point for such calculations [52].
Bearing this in mind, we derive bounds on the value of g e g N in the following way: Given experimental transition frequencies for several different intervals, e.g., ∆ exp b1a1 , ∆ exp b2a2 , ∆ exp b3a3 , etc., we calculate a value of R and a value of R p by matching these results with the corresponding theoretical frequency intervals, The values of these two parameters are determined by correlated χ 2 -fitting. We then obtain bounds on the coupling constant by finding the most positive and most negative values of g e g N for which the model fits the data at the 5% confidence level. The sensitivity to new physics arises because of the dependence of the NP shift on the quantum numbers n and l illustrated in Figs. 1-3. Put simply, states with high values of n and l are only weakly sensitive to new physics, whereas the opposite is the case for low-lying states.
Before describing the results of this analysis, we briefly discuss the existing experimental results relevant for this calculation and the related theoretical uncertainties. Further details about the calculation can be found in Appendix B 1.

A. Existing spectroscopic data for hydrogen
Clearly, detecting a NP interaction from spectroscopic data sets a challenging level of precision and accuracy on the measurements. Apart from the hyperfine splittings of the 1s and 2s states, which are not directly relevant here, the only hydrogen frequency intervals currently known to an accuracy better than 1 kHz are the 1s -2s interval, which has been measured with an experimental error of 10 Hz (i.e., a relative error of 0.004 ppt) [53,54], and intervals between circular states with n ranging from 27 to 30, for which unpublished measurements with an experimental error of a few Hz (about 10 ppt) have been made [55]. Circular states are states with |m| = l = n − 1.
The recommended value for the Rydberg constant is based on the 1s -2s measurement as well as on a number of measurements with a larger error [49,56]. The latter include measurements of the 2s -8s, 2s -8d and 2s -12d intervals made in the late 1990s with an experimental error ranging from 6 to 9 kHz (i.e., of the order of 10 ppt) [57][58][59]. Until recently, no other transitions between hydrogen states differing in n had been measured with an error of less than 10 kHz. However, the centroid of the 2s -4p interval has now been determined with an error of 2.3 kHz [60], and that of the 1s -3s interval with an error of 2.6 kHz (1 ppt) [61].

B. Theoretical uncertainty
The overall uncertainty on the SM predictions of hydrogen energy levels is mainly contributed by uncertainties on the values of the Rydberg constant, of the proton radius and of various QED corrections. Uncertainties on the values of other fundamental constants also contribute, although not in a significant way at the level of precision these energy levels can currently be calculated.
The uncertainty on the values of the Rydberg constant and the proton radius does not affect our calculation of the bounds on g e g N (recall that within our approach, these values are determined together with the bounds themselves in a self-consistent way).
Recent compilations of the relevant QED corrections and their uncertainties can be found in [50] and [51]. These corrections roughly scale as n −3 and strongly depend on l. Ref. [50] gives the combined theoretical uncertainty on the energy of a state of principal quantum number n as (2.3/n 3 ) kHz for l = 0, excluding the error contributed by the uncertainty on R p , and as less than 0.1 kHz for l > 0. Except for the 1s -2s interval, the experimental uncertainty rather than the theoretical uncertainty is thus the main limitation for setting bounds on g e g N based on the current spectroscopic data.

C. Bounds based on existing data
Bounds on the NP interaction strength derived as explained above are presented in Figs. 4(a) and 4(b), respectively for attractive and repulsive interactions. These results are based on three different sets of data, which we refer to as sets A, B and C. Set A groups all the existing high precision spectroscopic measurements in hydrogen, namely all the 18 experimental hydrogen transition frequencies included in the CODATA 2014 least square fit [49], the recent results of Ref. [60] for the 2s -4p interval and of Ref. [61] for the 1s -3s interval, and the results of Ref. [55] for the transitions between high lying circular states. The other two sets are the same as  Set A but without the 2s -4p results (Set B) or without the circular states results (Set C). The corresponding values of R obtained when assuming no NP shift are given in Table II, together with the recommended value of this constant [49], values based on the recent measurements of either the 2s -4p or the 1s -3s intervals [60,61] and a value based entirely on measurements of transitions between the circular states [55]. As is well known, the results of Ref. [60] are discrepant with both the CODATA results and those of Ref. [61] in regards to the values of R and R p , but yield a value of R p in good agreement with measurements in muonic hydrogen [62]. The values of R obtained from Dataset B are in close agreement with the CODATA 2014 value and have an uncertainty of a similar magnitude, although the CODATA fit also included spectroscopic measurements in deuterium and scattering data. Including the results of Ref. [60] in the fit reduces R significantly (the change is large because of the particularly small experimental error on these measurements).
Our main results for the current bounds on g e g N are based on Dataset A and are represented by the shaded areas in Fig. 4. They set a constraint of better than 10 −11 over the range of 10 1 -10 3 eV. As seen from the figure, the shape of the excluded area somewhat differs between attractive and repulsive interactions, particularly in the region around 100 eV. This difference indicates that the range of allowed values of g e g N is not centred on zerothough we emphasise that a value of zero remains compatible with the experimental data. The regions below the shaded areas indicate the range of values of g e g N compatible with the data, given the experimental and theoretical errors [63] Next we consider the effect of removing individual measurements from the calculation. Removing the recent 2s -4p measurement [60] has a considerable effect, not only weakening the overall bound, as expected, but also changing the shape of the excluded region. These differences reflect the aforementioned inconsistencies in the values of the Rydberg constant and the proton radius derived from the results of Ref. [60] with those obtained in the CODATA 2014 fit. The effect of removing this measurement illustrates the perils of selectively setting bounds using individual measurements or combinations of measurements. Whilst individual measurements may be precise, their accuracy can only be gauged against other measurements, particularly independent measurements of the same transitions.
Instead of removing the 2s -4p measurements, we now remove the unpublished circular state measurements of Ref. [55] and use Dataset C. The result is a substantial weakening of the NP bound for lower masses, illustrating the importance of using measurements of states with a large spatial extension when probing for a NP interaction with a low value of m X0 [8]. Although small, the NP shift of the circular states is not negligible when the range of the interaction is long enough. This leads to a decrease in the relative shift of these states compared to the low lying states when m X0 → 0, and hence to a weakening of the bounds on |g e g N | [64].
In summary, we have derived global NP bounds based on all available measurements for hydrogen, with no input from other atomic species. The sensitivity of the bound to individual measurements and to the Rydberg constant illustrates that bounds set using measurements on individual transitions should be treated with a degree of caution. The strong additional constraint provided by high-lying states at low masses motivate precision mea-surements for states with both higher n and l. For the latter we note the proposal of the Michigan group [65].
Before closing this section, we note that bounds on g e g N can also be found by comparing the values of R derived from different sets of transitions. For example, let R C (m X0 , g e g N ) and R D (m X0 , g e g N ) be the NPdependent values of R obtained by fitting the theoretical model respectively to Dataset C and to the circular state results of Ref. [55]. These two sets of data are completely independent of each other, and by contrast to R C (m X0 , g e g N ), the calculation of R D (m X0 , g e g N ) is insensitive to uncertainties on the proton radius and to poorly known QED corrections. The corresponding errors on these Rydberg frequencies, σ C and σ D , are also functions of m X0 and g e g N . As these errors are not correlated with each other, bounds on the NP coupling constant can be obtained by finding the most positive and most negative values of g e g N such that for a given choice of f (this constant sets the confidence limit of the bounds -we take f = 2). The results are also shown in Fig. 4 (the dotted curves). Cancellations of NP shifts are at the origin of the large weakening of these bounds between 1 and 10 keV. They are similar, below 300 eV, to those obtained from the global fit of the same set of data (the shaded areas). Compared to a global fit, however, this approach to setting bounds is potentially more sensitive to systematic errors in some of the measurements. We thus prefer to take the shaded areas as the best representation of the constraint on g e g N that can be set on the basis of the current body of spectroscopic work in hydrogen.

V. SCOPE FOR TIGHTER BOUNDS
Three factors limit the strength of the current bound shown in Fig. 4. The first is the experimental uncertainty of the measured energy levels. So far, only the 1s -2s interval has been measured with a relative uncertainty below the 0.01 ppt level. For higher states such as the measurements at n = 12, the ∼ 1 kHz uncertainty is approximately one hundred times larger or more. The second factor is the range of quantum numbers n and l for which precise data exist. The importance of additional measurements is highlighted in Fig. 4. Lastly, the limitations on the SM calculation of the energies also plays an important role. Here also there is much to be gained by working with higher-lying Rydberg states. In this section we consider the prospects for improvements in each of these three areas.

A. Improved measurements
In this section we consider the effect of reducing the current experimental uncertainty approximately 100fold, such that all transition frequencies in the dataset are known to the 10 Hz level currently available for the 1s -2s interval. As an aspirational goal we also consider what could be achieved with measurements at the 1 Hz level. A detailed discussion of future experiments is outside the scope of this article. Here we briefly discuss the dominant sources of uncertainties with the 10 Hz goal in mind. The focus is on laser spectroscopy of low-l states; improved measurements of circular Rydberg states are considered in [65].
The current measurement uncertainty includes contributions from both the background electromagnetic environment and atomic motion. Fundamental limits are provided by the radiative linewidth and black-body radiation (BBR). We calculated the radiative width and black-body shift and broadening of the relevant states (Appendix C). At n = 9, the radiative linewidth (which varies as n −3 ) is approximately 100 kHz for the s state and roughly ten times larger for the d state. The simple lineshape when radiative broadening dominates should enable line centres to be determined with high accuracy, with recent measurements in hydrogen determining line centres to one part in 10,000 of the linewidth [60]. As described in Appendix C, we find that black-body related uncertainties can be neglected even at 300 K provided that the temperature can be stabilised to 0.01 K.
Concerning stray fields, we note that the magnetic moment of low l states does not vary with n. Therefore, methods developed for precision measurements with low n states can be applied. For s-states the very small differential Zeeman shift is easily controlled at the sub-Hz level [59,66], while for d-states differential measurements such as those routinely carried out in optical atomic clocks [67] can be used to largely eliminate magnetic field errors. A much greater challenge is presented by the DC Stark shift, which scales as n 2 and n 7 for the linear and quadratic components respectively. A detailed analysis of the effect of the DC Stark shift on the hydrogen Rydberg spectrum is provided in [59]. In their experiments a stray field of ∼3 mV cm −1 was reported, leading to a final contribution to the uncertainty at the kHz level. However other experiment have shown that stray fields can be reduced to the 30 µV cm −1 level by performing electrometry with high-n states (n > 100) [68,69]. Drift rates as low as 2 µV cm −1 h −1 have also been measured [70]. Such measurements could be performed independently using a co-electrometry with a different species [69,70]. For a field of 30 µV cm −1 , the quadratic Stark effect is dominant for s-states, and measurements with 10 Hz uncertainty should be possible up to n = 23, with n ≈ 40 accessible if the stray field is determined to 1 µV cm −1 . For d-states, the linear Stark effect dominates, but differential measurements between different |m| states should enable the first order shift to be cancelled. The result-ing uncertainty thus becomes dominated by the residual quadratic shift.
Considering motional effects, we note that all measurements of hydrogen energy levels to date have been performed in atomic beams, where second-order Doppler effects limit the achievable linewidth to approximately 1 MHz. A complex velocity-dependent lineshape analysis is thus required to extract the true line center to the current 1 kHz accuracy [59].
In other atomic species, using ultracold atoms has enabled a dramatic reduction in the uncertainty of optical frequency measurements. Sub-10 Hz uncertainty has been achieved with untrapped atoms [71], while measurements based on atoms confined in magic-wavelength traps are entirely limited by the uncertainty in the microwave-based definition of the SI second [72,73].
For Rydberg states, experiments with ultracold atoms are dominated by the large level shifts due to the longrange van der Waals interaction [74], which scales as n 11 . Control over the number of atoms and interparticle distance and geometry is therefore essential. Confining atoms to a volume of ∼ 1 µm 3 would also largely eliminate errors due to field gradients. Therefore, a suitable platform could consist of individual hydrogen atoms confined in a single optical tweezer or tweezer array. Single-atom arrays have now been achieved with a growing range of atomic [75][76][77][78] and even molecular [79,80] species. Substantial hurdles exist for realising a similar system in hydrogen, not least the difficulty of laser cooling [81], which has so far proven essential for loading the optical tweezers. However alternative approaches such as loading from a hydrogen Bose-Einstein condensate [82], careful dissociation of laser-cooled hydride molecules [83] or in-trap Sisyphus cooling [84] may also provide possible routes. Here we assume that such a system may be realised, and that the contribution of the Doppler and recoil effects can be reduced below the natural linewidth of the transition by using well established two-photon spectroscopy techniques [59], possibly in combination with resolved sideband cooling [85]. Trap-induced AC Stark shifts are eliminated by extinguishing the trap light during the spectroscopy, as is common in Rydberg experiments with tweezer arrays.
Overall, we consider that a target of extending the range of states measured with an absolute uncertainty of 10 Hz or better to the full Rydberg series of s-and d-states up to a principal quantum number of n ≈ 40 is feasible. We note that this is still some way off the spectroscopic state-of-the-art achieved with cold trapped atoms. For circular states, 10 Hz uncertainty has already been achieved [55]; here achieving a precision of 0.1 Hz in future measurements seems feasible.

B. Improved theory
Improved measurements at the 10 Hz level would also provide a challenge to the current theory of SM correc-tions to hydrogen energy levels. Uncertainties in R and R p could be removed by using the global fitting procedure described in Section IV. Concerning the remaining correction due to QED and other effects, we note that the current uncertainty on the Lamb shift of the 2p 1/2 state is 21 Hz, including the uncertainty on the shift of the centroid of that level due to the hyperfine coupling [51]. As the theoretical error on QED and hyperfine corrections scales roughly like 1/n 3 and has been found to be smaller for states with larger orbital angular momentum, the theoretical error for the states with l > 0 is already expected to be below 10 Hz for n ≥ 3 and below 1 Hz for n ≥ 6. The situation for s-states is less clear. Current work assumes that the error on these corrections scales as n −3 , at least down to the 100 Hz level [50,51]. Given that the current theoretical uncertainty on the energy of the 2s state is about 2 kHz, achieving an accuracy of 10 Hz would require the evaluation of QED corrections that are currently rather poorly known. Alternatively, the data may be fitted to a theoretical model which does not rely on accurate values of the Lamb shift but instead treats the theoretical error on this quantity as a fitting parameter, assuming a n −3 scaling. We used such a model to obtain the illustrative results presented in Section V C (the method is outlined in Appendix B 2). However, further theoretical work would be necessary to confirm that the n −3 scaling still holds down to errors as small as 10 Hz or less. Each of the bounds shown in Fig. 5(a) was obtained by comparing the predictions of the Standard Model to a set of hypothetical data, the latter having been generated from a model including a NP shift. The details of the calculation are given in Appendix B 2.

C. Numerical illustration
The two blue curves plotted in Fig. 5(a) represent the bounds derived in this way from an arbitrary and hypothetical set of eight transitions between s-states, namely the 1s -2s, 2s -5s, 2s -8s, 2s -9s, 2s -11s, 2s -15s, 2s -21s and 2s -30s transitions. As seen from the figure, these results would improve the current spectroscopic bounds by two orders of magnitude over a wide range of values of m X0 , assuming an experimental error of 10 Hz. Reducing the error to 1 Hz would yield a three orders of magnitude improvement.
Using only transitions between states with l > 0 would remove the uncertainty on how the theoretical error scales with n. In practice, an experimental value for such a transition could be obtained, e.g., by measuring the 2s to (n, l) and 2s to (n , l ) intervals and subtracting one from the other to find the (n, l) to (n , l ) interval. The two orange curves plotted in Fig. 5(a)  the bounds derived from a set of transitions between dstates only (namely the 8d -9d, 8d -11d, 8d -15d, 8d -21d and 8d -30d transitions). While proceeding in this way has the advantage of avoiding the scaling issue, it has the disadvantage of taking into account only states with a relatively small NP shift. Correspondingly, and as is illustrated by the numerical results of Fig. 5(a), the bounds derived from such a set of data would be less stringent than those derived from data that include transitions from or between deeply bound states.
As mentioned above, bounds on g e g N can also be obtained by comparing the values of Rydberg constants derived from different sets of data. Assuming a 10 Hz experimental error and performing this comparison between the same sets of transitions as in Fig. 5(a) gives the bound represented by a solid curve in Fig. 5(b). This bound is slightly tighter but generally differs little from that obtained directly from the fit of the transitions between s-states. The dashed curve and dotted curve show that this bound could be lowered still fur-ther by also comparing these two values of the Rydberg constants with the value derived from transitions between circular states -i.e., transitions of the form (n, l = n − 1) ↔ (n = n + 1, l = n − 1). We consider two different sets of such transitions in Fig. 5(b). We took n = 10, 15, 20, 25 or 30 and assumed an experimental error of 0.5 kHz on these transitions to calculate the bound represented by a dashed curve, whereas for the bound represented by a dotted curve we took n = 40, 41, 42, 43 or 44 and assumed an experimental error of 0.1 Hz. Because the electronic density is concentrated further away from the nucleus when n > 40 than when n ≤ 30, adding the first or the second of these two sets of transitions lowers the bound in different ranges of values of m X0 .

VI. SUMMARY AND CONCLUSIONS
In summary, we have considered how the entire set of currently available spectroscopic data may be used to set global constraints on NP models that can be parameterized as a Yukawa-type interaction. Such interactions would naturally lead to so-called fifth forces which are a being searched for intensively [5,86].
Light force mediators have been intensively tested in lab experiments, e.g. through the Casimir effect [87]. As such searches rely in general on all atoms in a macroscopic object to contribute coherently and in concert to the resulting force on a test object, they do not probe directly the existence of a force on a microscopic level. This leaves large classes of new physics models untested. For example, forces mediated via kinetic mixing between a photon and a new Z [32,33] can easily avoid such bounds, as the atom as a whole is not charged under the fifth force. The experiments discussed above, however, would remain sensitive to such an interaction.
In addition, while other laboratory based experiments lose sensitivity for mediator masses above 100 eV, atomic spectroscopy for hydrogen atoms retains a good sensitivity up to masses of 10 keV. Thus, to our knowledge, the presented predicted limits provide the strongest constraints in laboratory based experiments obtained so far for that mass range.
The bounds we obtained in this work appear to be weaker than those set by astrophysical bounds. However, astrophysical bounds rely on the thermal production of light force mediators in stars [10,88]. Particles like chameleons avoid such production and thereby constraints from measurements of the energy transport in stars. Here atomic spectroscopy can help to close gaps in the landscape of Standard Model extensions and provide an independent test of the physics models underlying the assumptions of the models for the evolution of stars.
We further argue that this type of laboratory-based bound is unique, since it is independent of any manybody physics effects, such as astrophysical models or the complex subtleties of isotope shifts in many-electron atoms. Global constraints of this type also reduce the sensitivity to systematic errors in individual measurements, such as those which are currently giving rise to the so-called proton radius puzzle.
We therefore argue that there is a strong case for improved measurements in hydrogen based on extensions of current methods for precision optical frequency measurements in laser-cooled and trapped atoms. An important element would be extending the reach of measurements to higher principal quantum numbers, which has substantial benefits due to the dependence of the new physics shift on the shape of the wave function discussed in Section III. The ideal platform would be trapped single atoms or arrays of atoms with a well-controlled spacing, such as an optical tweezer array, opening also the tantalising prospect of an engineered many-body quantum system with a complete SM description.
An extension of this work would consider other simple atoms which have a complete SM description, such as D and He, or even positronium, muonic hydrogen or muonium. More sophisticated statistical analysis methods might enable measurements in all of these systems to be combined into a highly robust extended extended bound, or to create sensitive differential searches. Concrete limits could be obtained for various classes of new physics models, e.g., chameleons or kinetic mixing. Eq.
(2) being written in natural units, m X0 is expressed as an energy and r as the inverse of an energy. In atomic units, we have, instead, where the distance r is expressed in units of the Bohr radius a 0 , C in units of a −1 0 and B in units of α 2 m e c 2 a 0 , with m e the mass of the electron and α the fine structure constant. If r and r refer to the same point of space and r is expressed in eV −1 , then r = r c with the product c expressed in units of eV a 0 (the product c has the physical dimensions of an energy times a length and has a numerical value of 1 in natural units). Moreover, as m X0 r ≡ Cr if r and r refer to the same position, we see that m X0 and C are related by the equation That is, To relate the constant B to g e g N /4π, we note that in natural units g e g N is a pure number. However, since V (r) is actually an energy and r a length, Eq. (2) should really be written as Thus B, in atomic units, is (g e g N /4π) c with c expressed as a multiple of the product E h a 0 . Since E h a 0 = α c, c = (1/α)E h a 0 . This gives, to 10 s.f., Appendix B: Details of the fitting procedure

Bounds derived from current data
The calculation is outlined in Section IV. The full experimental dataset (set A) consists of the 18 measurements labelled A26.1 to A40.2 in Ref. [56], supplemented by the results of Refs. [60,61] and by a transition frequency for the transition between the n = 27 and n = 28 circular states calculated from the value of R quoted in Ref. [55]. (This value of R was derived from a small set of measurements of the n = 27 to n = 28 and n = 29 to n = 30 transitions, the former weighting more in the determination of R than the latter. No recommended value for either of these two transition frequencies is given in Ref. [55].) The set of data includes a measurement of the 2s 1/2 -2p 1/2 Lamb shift [89] recently reanalyzed in [90]. We use the revised value of this experimental result rather than its original value.
The correlation coefficients between the 18 measurements mentioned in Ref. [56] are given in that reference. We take the errors on these measurements to be uncorrelated with the errors on the measurements of Refs. [55,60,61] and the latter to be uncorrelated with each other.
The calculation of the terms R 2 p∆ ns biai and ∆ oc biai includes all the QED and hyperfine corrections listed in [50], as given in that paper and with additional input from [49,[91][92][93][94]. Following [50], we set the nuclear recoil contribution ∆ RR equal to 0 ± 10; taking ∆ RR = π/3 as proposed in [95] does not affect the results at the 0.1 kHz level of precision we are working here. As in [50], we assume a theoretical error of (2.3 kHz) h/n 3 on the s-state energies and of zero on the energies of the other states. The theoretical errors on the terms ∆ oc biai are thus completely correlated. The theoretical error on the factors ∆ g biai and∆ ns biai is too small to be relevant in the present context.

Projected bounds
We assume that each of the measured transition frequencies included in the set can be written as a sum of the form within experimental error, where R 0H = R 0 (m r /m e ) with R 0 the true value of the Rydberg frequency and ∆ corr ba is the sum of all the QED, hyperfine and other corrections predicted by the Standard Model. We take R 0H and ∆ corr ba to be the exact values of these quantities. We equate each of the experimental intervals to its Standard Model prediction, R H (1/n 2 a −1/n 2 b )+∆ corr ba +α corr ba , where R H is an effective Rydberg frequency obtained by fitting theory to experiment and α corr ba represents the theoretical error on ∆ corr ba . (This last term thus accounts for the error introduced by the uncertainty on the value of R p as well as the errors on the values of the QED and other corrections not calculated to a sufficient precision. We do not need to consider the uncertainty on the mass ratio (m r /m e ) separately from the uncertainty on R since this ratio is subsumed into the fitting parameter R H .) We assume that the n −3 scaling mentioned in Section V B holds for the s-states of interest, and that the theoretical error on the states with l > 0 is negligible Doing so for each of the N transitions of a same set yields the following overdetermined system: where A is a constant, and α exp biai represents the experimental error on the corresponding transition frequency. Simplifying these equations gives with δR H = R H − R 0H . We determine the NP bounds by finding the range of values of g e g N within which the left-hand sides of these equations fit the right-hand sides at the 5% confidence level, treating δR H and A as fitting parameters. As the experimental errors on different transitions measured using a same methodology could be expected to be mildly correlated, We assume a correlation coefficient of 0.1 between the experimental errors on the measured transition frequencies belonging to a same set of data. The bounds shown in Fig, 5(b) follow from Eq. (14), with the Rydberg frequencies R C and R D replaced by the corresponding values of δR H . For simplicity, we assumed no correlation between the errors on these quantities.

Appendix C: Black body radiation
We have calculated the BBR shift of the Rydberg states of interest in order to ascertain the precision on the thermometric measurements required in our approach. Our results complete those of Refs. [96] and [97], which do not extend high enough in principal quantum numbers. Except where specified otherwise, we use atomic units throughout this appendix.
To second order in the electric field component of the BBR field, the thermal shift of a state a at a temperature T can be written as [96][97][98] where the r i 's are three orthogonal components of the electron's position operator, k is Boltzmann constant, the summation over b runs over all the atomic states dipolecoupled to the state a, ω ab = E a − E b where E a and E b are the energies of the respective states, and The BBR field also depopulates state a by inducing transitions to other states at a rate approximately equal to Γ BB a , where with U (y) = |y| 3 /(exp |y| − 1) [96][97][98]. Γ BB a does not include losses due to the BBR-induced Stark mixing of degenerate states of opposite parity, which is significant in hydrogen [97]. Eqs. (C1) and (C3) also neglect nondipolar transitions and corrections of fourth order in the BBR electric field [99][100][101]; their contributions are small and can be neglected for our purposes. Local anisotropies of the BBR field may also need to be factored in when comparing to experiment [102].
We evaluate F (y) by contour integration in the complex x-plane. This method bypasses the need of a careful treatment of the singularity at x = ±y inherent in the direct calculation of a Cauchy principal value and only involves straightforward numerical quadratures. Namely, we introduce the complex function where the integration contour starts at z = 0 and goes to Re z → ∞ in the lower half plane, avoiding the zeroes of exp(z) − 1. In practice, we use a rectangular contour running from z = 0 to z = −i/2 and from z = −i/2 to z = 50 − i/2, which is well adapted to the range of values of y involved in this work. We have F(y) ≡ Re F C (y) owing to the relation and moreover We calculate the dipole matrix elements b|r i |a by solving Eq. (8) for each of the states a and b. Having the corresponding generalized eigenvectors, c a and c b , we obtain b|r i |a as c † b R i c a , where R i is the matrix of elements S * n l b (r)Y * l b m b (θ, φ) r i S n la (r)Y lama (θ, φ) d 3 r.
Substituting these results into Eq. (C1) gives the shift in the non-relativistic approximation. We correct this for spin-orbit coupling by replacing the non-relativistic angular factors by the appropriate expressions [96] and evaluating the Bohr transition frequencies ω ba using the relativistic energies. The summation over the intermediate states b runs over all the generalized eigenvectors of the matrix H 0 of the relevant symmetry, including those corresponding to positive eigenenergies. Doing so ensures (assuming that the basis is large enough) that the shifts and widths properly include the contribution of the continuum, which can be significant [103,104]. We use 300 Sturmian functions for each symmetry. While the BBR shift depends to some extent on the hyperfine structure of the levels [105], taking it into account would not affect the results at the level of precision required by the present investigation. The calculations of radiative widths mentioned in the text use exactly the same numerical method in regards to the computation of the required dipole matrix elements.
The BBR shift of hydrogen state may be significant compared to the NP shift at the relevant values of g e g N , and may even be considerably larger. For example, for n = 10 and l = 0, |δ NP nl | is at most 0.7 kHz when g e g N = 1 × 10 −12 whereas δ BB nl is approximately 1.1 kHz at 300 K [96]. At least in principle, this shift can be removed from spectroscopic data for hydrogen since it can be accurately calculated for this atom. In practice, however, taking it correctly into account requires a sufficiently precise determination of the temperature of the BBR field at the location of the atoms, and perhaps also of its inhomogeneity and its deviation from of an ideal Planck distribution. In situ temperature measurements with an uncertainty of the order of 0.01 K have been achieved using platinum resistance thermometers [106]. Spectroscopic measurements of Rydberg states have also been proposed to determine the temperature of the BBR background with a similar uncertainty [107]. The BBR energy shift of the high Rydberg states is approximately π(kT ) 2 /3c 3 [96,108]. An error of 0.01 K on T at 300 K translates into an error of 0.2 Hz (or lower) on the BBR shift of these states. However, this error on the temperature would have a smaller impact on measurements of the energy difference between Rydberg states made in a same apparatus because they all shift by roughly the same amount. This point is illustrated by Fig. 6, which shows the accuracy to which T must be known to reduce the error on the BBR shift to less than 0.2 kHz in measurements of transitions between sstates made at room temperature. An easily achievable accuracy of 0.5 K is sufficient for transitions between the lowest states or between high Rydberg states (the former because they shift little, the latter because they shift similarly). The requirements are more stringent for transitions between relatively low lying states and high Rydberg states, particularly for low lying states with n ≈ 5 (whose shift is larger and of opposite sign to that of states with n 5 [96]). It should be noted that the error that can be tolerated on T roughly scales with the maximum error on the transition frequencies. If an accuracy of 10 Hz is sought, rather than 0.2 Hz, knowing the temperature would not need to be known to better than 0.5 K for any frequency interval.