Probing QED and fundamental constants through laser spectroscopy of vibrational transitions in HD+

The simplest molecules in nature, molecular hydrogen ions in the form of H2+ and HD+, provide an important benchmark system for tests of quantum electrodynamics in complex forms of matter. Here, we report on such a test based on a frequency measurement of a vibrational overtone transition in HD+ by laser spectroscopy. We find that the theoretical and experimental frequencies are equal to within 0.6(1.1) parts per billion, which represents the most stringent test of molecular theory so far. Our measurement not only confirms the validity of high-order quantum electrodynamics in molecules, but also enables the long predicted determination of the proton-to-electron mass ratio from a molecular system, as well as improved constraints on hypothetical fifth forces and compactified higher dimensions at the molecular scale. With the perspective of comparisons between theory and experiment at the 0.01 part-per-billion level, our work demonstrates the potential of molecular hydrogen ions as a probe of fundamental physical constants and laws.


Introduction
The possibility of accurate quantum electrodynamics (QED) calculations and the presence of narrow optical transitions between long-lived vibrational levels make laser spectroscopy of H 2 + and HD + a sensitive tool to test fundamental physics at the molecular scale.Examples include tests of relativistic quantum mechanics and QED 1,2,3,4 , and searches for physics beyond the Standard Model at the molecular scale 5,6 and beyond General Relativity 7,8 .
Furthermore, spectroscopy of the molecular hydrogen ion has long been proposed as a means to determine the value of fundamental physical constants 1,2,4,9 , and to realise precise molecular clocks 8,10 .
Nearly all the above applications involve a comparison of state-of-the-art molecular theory with accurate results from (most often laser) spectroscopy.The most accurate results so far have been obtained using HD + ions stored in radiofrequency (r.f.) traps, sympathetically cooled by laser-cooled Be + ions 3,11 , allowing tests at a relative uncertainty of a few parts per billion (p.p.b.)However, the most precise experimental result 3 , a frequency measurement of the (v,L): (0,0) -(1,1) rovibrational line at 5.1 m in HD + with a relative uncertainty of 1.1 p.p.b., was found to disagree by 2.7 p.p.b. with a more accurate theoretical value obtained from state-of-the-art ab initio molecular theory 1,2 .This disagreement has thus far been unresolved, and additional high-precision experimental data are needed to draw conclusions about the validity of the theoretical framework, and to open up the wide range of applications of molecular hydrogen ion spectroscopy mentioned above.
Here, we report on an optical frequency measurement of the (v,L): (0,2) -(8,3) vibrational overtone transition in HD + at 782 nm.Our experimental result is found to be in agreement with the theoretical prediction to within 0.6(1.1)p.p.b., thereby confirming the validity of relativistic quantum mechanics and QED in a vibrating molecular system at an unprecedented level.We subsequently exploit the agreement between theory and experiment for the first determination of the proton-electron mass ratio from a molecular system, and to put tighter constraints on the strength and range of 'fifth forces' at the molecular scale.

Results
Experimental procedure.Our experimental apparatus and procedure for HD + spectroscopy are described in detail in the Methods section.In brief, trapped HD + molecular ions are cooled to ~10 mK by storing them together with laser-cooled Be + ions, and we excite and detect the (v,L): (0,2) - (8,3) transition by resonance-enhanced (1+1') multiphoton dissociation (REMPD) spectroscopy, see Fig. 1a.We acquire a spectrum by observing the loss of HD + due to REMPD, inferred from the change in the monitored Be + fluorescence induced by r.f.excitation of the HD + secular motion, for different values of the 782 nm spectroscopy laser frequency (Methods and Fig. 2a).The latter is stabilized and counted using an optical frequency comb laser (Methods).
Line shape model for nonlinear least-squares fitting.Hyperfine interactions lead to a manifold of 83 lines (Figs. 1b,c), the ~25 strongest of which are located within the range (-110 MHz, 140 MHz) around the hyperfine-less rovibrational frequency,  0 .Each hyperfine component is Doppler broadened to ~16 MHz, and the hyperfine structure is only partially resolved.Using a realistic spectral line shape function (Methods), we employ nonlinear least squares fitting to extract not only the transition frequency  0,fit , but also the intensity of the 782 nm laser, I L , the motional temperature of the HD + ions, T HD+ , and the temperature of the Be + ions during secular excitation, T i , which are relevant parameters for the spectral analysis (see Methods).The 782 nm intensity we find agrees (within the fit error) with the intensity estimated from the laser beam waist and beam alignment uncertainty.Likewise, the fitted temperature agrees well with the results of molecular dynamics (MD) simulations.
The relevance of the additional fit parameters is illustrated by the correlation matrix, , of the estimated parameters, which reveals significant correlations between  0,fit and T HD+ , and  0,fit and I L (Fig. 2c).Adding fit parameters results in an increased error in the fitted value  0,fit , which rises from 0.25 MHz (0.65 p.p.b. relative to  0,fit ) to 0.33 MHz (0.85 p.p.b.) when T i , T HD+ , and I L are added as free fit parameters.Even so, the spectral fit result presented here marks the first time that a sub-p.p.b. resolution is achieved in molecular-ion spectroscopy.
Systematic effects and frequency value.The fit result  0,fit is corrected for various systematic frequency offsets due to electric and magnetic fields.We calculate the a.c. and d.c.Stark effect ab initio with high accuracy following the approach of Karr 10 to find the frequency shift due blackbody radiation (BBR), the r.f.trap field, and the electric fields of the lasers.The total Stark shift of -1.3(1) kHz is dominated by the shifts due to the 313 nm, 532 nm and 782 nm lasers (which are all on during excitation), with individual contributions of 0.008(1) kHz, -0.45(7) kHz, and -0.87(13) kHz, respectively.Here, the uncertainties of the laser beam intensities are the dominant source of frequency uncertainty.We also evaluate the Zeeman effect, which stems from the static quantization field of 0.19 mT directed along the z-axis of the trap (Methods).The polarization of the 782 nm laser is such that it induces   and   transitions at equal rates.The Zeeman effect leads to a shift to  0,fit by -16.9(3.2) kHz, the uncertainty of which is due to a possible 2% imbalance between the   and   rates caused by imperfect polarization optics and depolarization due to the windows of the vacuum chamber.Another uncertainty stems from the accuracy of the theoretical hyperfine structure 12 , which enters through our spectral line shape function.We repeated the fit procedure with a spectral line shape function based on the hyperfine structure obtained with slightly altered spin coefficients (within their uncertainty range), for which we observe essentially no shift of  0,fit .Compared to the 0.33 MHz statistical fit uncertainty of  0,fit , the contribution of the above line shifts to the total frequency value and uncertainty is negligible, as are the frequency shifts due to the second-order Doppler effect 13 (<5 Hz) and the electric-quadrupole shift 14  (with ionic product kinetic energies of 0.25, 0.36, and 0.66 eV, respectively) into account here.We determine reaction rates from measured loss rates of Be + and HD + , which are in agreement with expected Langevin reaction rates given the pressure of 110 -8 Pa in our vacuum apparatus 16 .We include these processes in realistic MD simulations (Methods) employing leapfrog integration with an adaptive time step, to ensure that collisions between energetic ions are correctly handled.Under these conditions, our MD simulations reveal average HD + z-velocity distributions which deviate significantly from thermal (Gaussian) distributions, an effect which has hitherto been neglected by the widespread assumption that laser-cooled ion ensembles exhibit thermal velocity distributions.As shown in Fig. 3, we empirically find that the velocity distributions are better represented by q-Gaussians 17 , which have the additional advantage that the effect of a given chemical reaction can be parameterized by a corresponding q-value (with q=1 corresponding to a Gaussian distribution).Taking all reactions and reaction rates into account, we find that the velocity distribution is best characterized by a q-Gaussian with q ranging between 1.00 and 1.07, depending on the assumed thermalisation rate and the REMPD rate.We therefore include q-Gaussians in our spectral line shape function (Methods).Comparing scenarios with minimum and maximum expected thermalisation rates (corresponding to q=1.07 and q=1.00, respectively) we find a mean shift of -0.25 (25) MHz with respect to the case of q=1.00.The origin of this shift lies in the overlap and saturation of Doppler-broadened hyperfine components in the spectrum (Fig. 2d).A similar frequency shift may occur when micromotion sidebands are present.Here we make a distinction between excess radial micromotion caused by radial static offset fields, and axial micromotion which could arise from an axial r.f.field due to geometrical imperfections of our ion trap.A finite-element analysis of the trap's electric field reveals that such an axial field will be approximately constant over the extent of the Coulomb crystal, which implies that the corresponding micromotion contribution to the line shape will be homogeneous.We use the fluorescence correlation method of Berkeland et al. 13 to find a small micromotion amplitude of 11(4) nm amplitude along the 782 nm laser beam as follows.Measurements of the r.f.field amplitude are performed on the Be + ions using the 313 nm laser beam (which co-propagates with the 782 nm laser beam).Repeating such measurements for various values of the trap r.f.voltage allows separating the axial and the radial component (the latter being due to the residual projection of the laser beam onto the radial direction).We find that the axial micromotion component is dominant, which is explained as follows.Firstly, the 782 nm laser is aligned parallel to the trap axis, so that the radial micromotion amplitude is suppressed by the nearzero angle between the 782 nm laser and the radial directions.Secondly, since the surrounding Be + ions shield the HD + from static radial offset fields, the HD + ions are trapped symmetrically about (and close to) the nodal line of the r.f.field.These two conditions limit the contribution of the radial micromotion to well below the measured axial amplitude of 11(4) nm.The resulting sidebands are included in the spectral line shape function, leading to an additional frequency shift of -55 (20) kHz with respect to the case of zero micromotion.

Discussion
Our frequency value  0 of the (v,L): (0,2) -(8,3) transition is in good agreement with the more accurate theoretical value 2  0,theo = 383,407,177.150(15) MHz.The contribution of the QED terms 1,2,18 to this frequency amounts to -1547.225(15)MHz (-4035.46(4)p.p.b.), and our measurement therefore confirms the validity of QED in a molecular system at an unprecedented level of 2.7×10 -4 .We furthermore note that our measurement is sensitive to (and in agreement with) QED terms up to order m e  7 , given that the m e  7 term contributes 780( 15) kHz (i.e.2.03(4) p.p.b.) to the transition frequency.
Salumbides et al. 5 showed that spectroscopic tests of molecular QED can be used to set upper bounds on a hypothetical 'fifth force', acting between hadrons at separations of the order of 1 Å, and arising from the exchange of unknown, massive virtual particles.Modelling such a hadron-hadron interaction with a Yukawa-type potential of the form ħc   e r/ /r and computing the resulting frequency shift to  0 , we can exploit the 1 p.p.b. agreement between theory and experiment to rule out (at the 90% confidence level) interactions at a range  = 1-2 Å (corresponding to force-carrying particles with mass m Y = 1-2 keVc -2 ) with an interaction strength relative to the fine-structure constant,   /,larger than 5-8×10 -10 (Fig. 4).In a similar way, applying the Arkani-Hamed-Dimopoulos-Dvali formulation to probe the effect of compactified higher dimensions on energy levels in molecules 6 , we place an improved upper bound of 0.5 m on the compactification radius of higher dimensions in eleven-dimensional M-theory.
Four decades ago, Wing et al. 4 suggested that molecular theory could one day be used to is still 31 times higher 22 , we point out that our method yields  as a single parameter from molecular vibrations, whereas most other  values are ratios of individual determinations of the electron and proton relative atomic masses (exceptions are  determinations from atomic laser spectroscopy 23,24 ).Therefore, the agreement of  HD+ with all other values of  forms an additional consistency check of the various (and vastly different) methods used (Fig. 5b).In particular, it implies that relativistic quantum mechanics and QED consistently describe at the few-p.p.b. level such diverse systems as the bound electron 22,25,26 , antiprotonic helium 23 , and the molecular hydrogen ion.Furthermore, of all the methods which produce  as a single parameter, our approach is surpassed only by spectroscopy of antiprotonic helium , which is 2.3 times more precise but additionally requires charge, parity and time reversal invariance 23 (Fig. 5a).
In principle, the transition frequency   depends on the value of other fundamental constants as well, such as the deuteron-proton mass ratio 7 , m d /m p , the fine structure constant, and the proton charge radius 1,2 .However, the sensitivity of   to changes in  is known to be three times larger than the second largest one, the m d /m p sensitivity 7 .
Moreover, if we propagate the uncertainties of the 2010 CODATA values of the fundamental constants through the sensitivity relations, we find that the relative contributions by , m d /m p ,  and the proton radius to the theoretical uncertainty of  0 are 154, 11.6, 0.004, and 5.13 parts per trillion, respectively.We therefore conclude that  is the correct parameter to constrain.
The error budget in Table 1 reveals that the experimental uncertainty is limited by Doppler broadening.To overcome this, more involved Doppler-free two-photon spectroscopy of HD + has been proposed 9 , which should reduce the uncertainty to below 1×10 -12 .This should enable a comparison between theory and experiment which would not only test the QED description of chemically bonded particles at an unprecedented level, but also produce a new value of  with an uncertainty below 0.1 p.p.b., surpassing the most precise determination of  obtained from independent electron and proton relative atomic mass measurements 22 , which represent a completely different method (Fig. 5c).Our work demonstrates the potential of molecular hydrogen ions for the determination of mass ratios of fundamental particles, as well as stringent tests of QED, and searches for new physics.

Methods
Experimental procedure.We typically store 40 to 85 HD + ions in a linear rf trap, together with about 750 Be + ions, which are laser cooled to a temperature of 5-10 mK using a continuous-wave (CW) 313 nm laser beam propagating along the symmetry (z) axis of the trap.As a consequence, only the axial motion of the Be + ions is laser-cooled directly.
However, the three-dimensional extent of the Coulomb crystal ensures good coupling of the axial motion to the radial motion of the Be + and HD + ions, so that both Be + and HD + are efficiently cooled in all three dimensions.Although larger numbers of ions may be trapped, smaller ion numbers ensure better reproducibility of the experimental conditions.The stronger confinement of HD + by the pseudopotential leads to the formation of a string or zigzag structure of HD + ions along the nodal line of the r.f.field, which coincides with the z-axis.
Be + ions arrange themselves in a three-dimensional ellipsoidal structure surrounding the HD + .At temperatures below 0.1 K, the two-species ion ensemble solidifies into a Coulomb crystal.313 nm fluorescence photons emitted by Be + are imaged onto an electron multiplied charge-coupled device (EMCCD) camera and a photomultiplier tube.We obtain a measure of the number of trapped HD + ions by resonantly driving their radial secular motion (~800 kHz) using an a.c.electric field 19 .MD simulations indicate that this 'secular excitation' heats up and melts the ion ensemble, heating the Be + ions to an average temperature (in the zdirection) T Be+,av proportional to the number of HD + , N HD+ .Typically, T Be+,av = 2-4 K.For a fixed cooling-laser detuning =-18  (with =2×19.4MHz the natural linewidth of the 313 nm cooling transition), this temperature rise leads to Doppler broadening and, thus, to a significant rise in the 313 nm fluorescence rate with average value F. Whereas previous work approximated the fluorescence rise versus N HD+ by a linear dependence 27 , we realistically model the (nonlinear) dependence of F on N HD+ , and take this into account in our analysis (see Spectral line shape function for fitting).
To excite the (v,L): (0,2) -(8,3) rovibrational transition at 782 nm, we send a CW 782 nm laser beam along the z-axis, counter-propagating the 313 nm laser beam.The 782 nm radiation is obtained from a titanium:sapphire laser with a linewidth of 0.5 MHz, which is frequency locked to an optical frequency comb laser as follows.The frequency of the optical beat note of the 782 nm laser with a nearby mode of the comb is measured by a counter every 30 ms.From the measured beat-note frequency, the comb repetition rate, comb carrier-envelope offset frequency, and the comb mode number we determine the laser's optical frequency.The counter output is fed into a digital feedback loop, which controls the 782 nm laser so as to stabilise the beat-note frequency to a pre-set value.The comb is fully referenced to a GPS-disciplined rubidium atomic clock (providing long-term relative accuracy at the level of 2×10 -12 ), and the frequency uncertainty of our optical frequency measurement system (10 kHz) is limited solely by the frequency instability of the locked 782 nm laser averaged over the 10 s of REMPD.To detect excitation to the v=8 state by the 782 nm laser, we overlap this beam with a CW 532 nm laser beam, leading to REMPD of HD + (Fig. 1a).An experimental cycle looks as follows.We first load a fresh sample of HD + and find a measure, F i , of the initial ion number, N i , by secular excitation.We subsequently lower the 313 nm cooling laser intensity and detuning to reduce the ion temperature to ~10 mK, and expose the ions to the 313 nm, 782 nm and 532 nm lasers for 10 s.Afterwards, we apply secular excitation to determine the fluorescence level F f corresponding to the remaining ion number, N f , and we define a signal S = (F i -F f )/F i .Repeating this cycle for different values of the 782 nm laser frequency, , we obtain a spectrum consisting of 1772 points (Fig. 2a).
Noise in the spectrum is primarily due to the stochastic nature of REMPD of small HD + ensembles with a mean rotational-hyperfine occupancy of the order of one ion per state.
Spectral line shape function for fitting.Hyperfine rate equations allow computing the evolution of a sample of HD + ions, with an initial thermal rotational distribution 28 corresponding to T =300(5) K, under the influence of REMPD lasers, BBR and losses due to chemical reactions with background-gas molecules (Figs.1b,d).Population in rotational states with L=0-5 is included, thus ignoring the 2.4% population in L=6 and higher.Accurate hyperfine line strengths at the magnetic sub-state level are obtained by extending the approach of Koelemeij 29 so as to include hyperfine structure and the Zeeman effect.The response to the 782 nm laser of each hyperfine level is modelled using a Doppler-broadened profile based on q-Gaussians.For the assessment of the Zeeman effect, the hyperfine rate equations are adapted to include Zeeman line shifts and broadening due to hyperfine line splitting (which remains much smaller than the 16 MHz Doppler width).The q parameter is determined from realistic MD simulations (which take into account the time-dependent trap potential and momentum changes by scattering of photons from the cooling laser) with an uncertainty limited by the minimum and maximum expected rates of heating events.Solving the rate equations allows predicting the fractional loss of HD + , , as a function of the 782 nm laser frequency relative to the hyperfine-less frequency,  0 , the intensity, I L , of the 782 nm laser, and the HD + motional temperature, T HD+ .We obtain a smooth, continuous spectral line shape function for fitting as follows.First, we compute  on a dense three-dimensional grid of values (, I L , T HD+ ), which we subsequently interpolate to obtain a four-dimensional function ( 0,fit , I L , T HD+ ).Using the shorthand notation ( 0,fit , I L , T HD+ ), and making use of the linear dependence of T Be+,av on N HD+ found from MD simulations, the average Be + temperatures during secular excitation of HD + before and after REMPD are given by T i  T Be+,av (N i ) and T f  T Be+,av ((1-)N i ), respectively.We consecutively insert these temperature values in the scattering rate function g, defined as which includes the beryllium mass, m Be+ , and the fixed 313 nm laser detuning, wave vector, k, and saturation parameter, s.The integration is performed over the distribution of Be + velocities v along k, i.e. v k  k•v/|k|.This function is used to model the fluorescence rise F observed in the experiment, and allows us to construct a five-dimensional fit function with g 0 the steady-state scattering rate before secular excitation.The function G realistically models the signal S (Fig. 2c), and is used for fitting.
We employ the above model also to estimate the uncertainties due to several systematic effects.For example, MD simulations predict a slight increase of T HD+ and q with increased REMPD rates in the scenario of inefficient thermalisation of fast D + .We can estimate the frequency shift due to this effect by making both T HD+ and q REMPD-rate dependent.This leads to an additional frequency uncertainty of 61 kHz, which is included in the error associated with Doppler effects due to chemistry (Table 1).From our model we furthermore deduce that ignoring states with L6 introduces an uncertainty of 32 kHz, while the possible 5 K error in the BBR temperature estimate (which takes into account day-to-day temperature variations and the possibility that the trap electrodes may be at a slightly higher temperature due to r.f.dissipation) translates to a 5 kHz effect.Values labelled with an asterisk correspond to frequency shifts which are accounted for by the spectral fit function, rather than being subtracted from  0,fit .Values labelled with a dagger are non-zero but negligibly small and therefore ignored here.
translate measured vibrational frequencies of HD + to a new value of the proton-electron mass ratio, .The high accuracy of our result and the good agreement with theory, which assumes the 2010 Committee on Data for Science and Technology (CODATA) recommended value  CODATA10 , now allow the first determination of  from molecular vibrations.Previously 19 we had derived the sensitivity relation / = -2.66  /  , which we employ to adjust  CODATA10 to a new value,  HD+ , such that the theoretical frequency matches our experimental value.We thus find  HD+ =1,836.1526695(53),with a relative uncertainty of 2.9 p.p.b. which approaches that of the values taken into account in the 2010 CODATA adjustment20 .For example, the value reported by Farnham et al. is only a factor of 1.3 more precise than our result21 .While the precision of a very recent determination bySturm et al.

Figure 1|
Figure 1| Partial level diagram, excitation scheme and interaction with blackbody

Figure 2|
Figure 2| Spectral line shape function and fit to REMPD spectrum.a, Measured

Figure
Figure 3| Non-thermal velocity distributions.When kinematic effects of (laser-induced) chemistry are included in MD simulations, the HD + z-velocity distribution exhibits non-

Figure
Figure 4| Constraint on fifth forces between hadrons at the Ångstrom scale.The high

Figure
Figure 5| Determination of the proton-electron mass ratio.a, Comparison of values of 