High-precision spectroscopy of the HD+ molecule at the 1-p.p.b. level

Recently we reported a high precision optical frequency measurement of the (v,L):(0,2)->(8,3) vibrational overtone transition in trapped deuterated molecular hydrogen (HD+) ions at 10 mK temperature. Achieving a resolution of 0.85 parts-per-billion (p.p.b.) we found the experimental value ($\nu_0= 383,407,177.38(41)$ MHz) to be in agreement with the value from molecular theory ($\nu_\text{th}=383,407,177.150(15)$ MHz) within 0.6(1.1) p.p.b. [Biesheuvel et al., Nat. Commun. 7, 10385 (2016)]. This enabled an improved test of molecular theory (including QED), new constraints on the size of possible effects due to 'new physics', and the first determination of the proton-electron mass ratio from a molecule. Here, we provide the details of the experimental procedure, spectral analysis, and the assessment of systematic frequency shifts. Our analysis focuses in particular on deviations of the HD+ velocity distribution from thermal (Gaussian) distributions under the influence of collisions with fast ions produced during (laser-induced) chemical reactions, as such deviations turn out to significantly shift the hyperfine-less vibrational frequency as inferred from the saturated and Doppler-broadened spectrum, which contains partly unresolved hyperfine structure.


Introduction
Because of their relative simplicity, molecular hydrogen ions such as H + 2 and HD + are benchmark systems for molecular theory [1,2], and suitable probes of fundamental physical models [3,4]. Already four decades ago, Wing et al. [5] performed measurements of rovibrational transitions in HD + and conjectured that such measurements could be used to test quantum electrodynamics Send offprint requests to: j.c.j.koelemeij@vu.nl Present address: Statistics Netherlands (CBS), P.O. Box 24500, 2490 HA The Hague, The Netherlands (QED) theory in molecules, and lead to an improved value of the proton-to-electron mass ratio, µ. Today, QED corrections to ro-vibrational transition frequencies in H + 2 and HD + have been calculated up to the order m e α 8 , with m e the electron mass and α the finestructure constant, leading to relative uncertainties of the order of 3-4 × 10 −11 [1,2]. The frequencies of these transitions are typically in the (near-) infrared with linewidths below 10 Hz, which makes them amenable to high-resolution laser spectroscopy. Optical clocks based on trapped atomic ions have shown that laser spectroscopy at very high accuracy (below one part in 10 17 ) is possible [6]. Recent theoretical studies point out that also for molecular hydrogen ions experimental uncertainties in the 10 −16 range should be possible [7,8].
Recent frequency measurements of rovibrational transitions in HD + were done with a reported uncertainty of 1-2 p.p.b. [9,10]. However, the most precise measurement of a rovibrational transition in HD + so far resulted in a 2.7 p.p.b. (2.4σ) discrepancy with more accurate theoretical data [10]. The question whether this offset is a statistical outlier, or caused by an experimental systematic effect or possible 'new physics' has remained unanswered. Therefore, additional experimental data on HD + are needed. In this article we present the result of a high-precision frequency measurement of the rovibrational transition (v,L):(0,2)→ (8,3) in the HD + molecule, and compare it with state-of-the-art molecular theory to find excellent agreement. From this comparison we subsequently derive new constraints on the size of effects due to possible new physics in HD + , and we determine a value of µ for the first time from a molecular system [11].
This article is organized as follows. In Sec. 2 we briefly review the theory of HD + relevant to this experiment. In Sec. 3 we describe a setup for spectroscopy of trapped HD + ions, sympathetically cooled by laser-cooled beryllium ions. The analysis of the spectroscopic data, sys-tematic frequency shifts and the results are discussed in Sec. 4, followed by the conclusions in Sec. 5.

Calculation of ro-vibrational frequency transitions in HD +
The calculation of ro-vibrational energies in quantum mechanical three-body systems is usually divided into a non-relativistic part, complemented by the relativistic and radiative parts, and finite nuclear size contributions. The non-relativistic part is calculated through solving the three-body Schrödinger equation, which can be done with practically infinite precision [12] (up to a relative precision of ∼ 10 −32 [13,14]). The resulting wavefunctions allow an analytical evaluation of the Breit-Pauli Hamiltonian and the leading-order radiative corrections.
Recently, the precision of relativistic and radiative energy corrections to the non-relativistic energies was strongly improved. With the inclusion of the full set of contributions of order mα 7 and leading-order terms of order mα 8 , the relative uncertainty is now below 4 × 10 −11 . For example, the theoretically determined value of the HD + (v,L):(0,2)→(8,3) transition frequency, ν th , is 383,407,177.150 (15) MHz, and has a relative uncertainty of 4 × 10 −11 [1,2]. Note that the specified error (within parentheses) does not include the uncertainty of the fundamental constants used. By far the largest contribution is due to the 4.1×10 −10 uncertainty of the 2010 Committee on Data for Science and Technology (CO-DATA) value of µ [15], which translates to a frequency uncertainty of 59 kHz.

Hyperfine structure and rotational states
Since the HD + constituents possess nonzero spin, the ro-vibrational transition spectra contain hyperfine structure due to spin-spin and spin-orbit couplings. In [16] the hyperfine energy levels and eigenstates in HD + are calculated by diagonalization of the effective spin Hamiltonian H eff = E 1 (L · s e ) + E 2 (L · I p ) + E 3 (L · I d ) + E 4 (I p · s e ) + E 5 (I d · s e ) + E 6 K d (L, I p , s e ) + E 7 K d (L, I d , s e ) where the spin coefficients, E i , are obtained by averaging the Breit-Pauli Hamiltonian over the nonrelativistic wavefunctions, L is the rotational angular momentum operator, and s e , I p and I d are the electron, proton and deuteron spin operators. K d and K Q are spherical tensors composed of angular momenta, whose explicit form is given in [16]. The strongest coupling is the protonelectron spin-spin interaction (the term in E 4 in Eq. (1)), and the preferred angular momentum coupling scheme is F = I p + s e S = F + I d J = L + S. (2) For states with L = 0, L = 1 and L ≥ 2 this leads to the presence of 4, 10 and 12 hyperfine levels, respectively, which are distributed over an energy range of ∼1 GHz ( Fig. 2(b)). Diagonalization produces eigenstates φFS JM J , with the magnetic quantum number M J corresponding to the projection of J onto the laboratoryfixed z-axis. Note that after diagonalization the quantum numbersF ,S are only approximately good, and a hyperfine eigenstate can be expressed in the 'pure' basis states |F SJM J in the form with real-valued coefficients β F SJ . In [16] the hyperfine levels in vibrationally excited states in HD + are calculated with an error margin of ∼ 50 kHz. This uncertainty might propagate through a spectral analysis employing a lineshape model which includes the theoretical hyperfine structure (as we do below). However, the uncertainty contribution to the determined hyperfine-less (v,L):(0,2)→(8,3) transition frequency is expected to be at least five times smaller due to strong correlation and partial cancelation of the error (see also Sec. 4.6.3).

Determination of transition rates
Because of the hyperfine structure, the spectrum of the (0,2)→(8,3) electric dipole transition consists of a large number of hyperfine components. Together with Doppler broadening, this leads to an irregular lineshape. In addition, the excitation laser may address multiple hyperfine states simultaneously, which are furthermore coupled to other rotational states by the ambient 300 K blackbody radiation (BBR) field. Therefore, for the analysis of the (0,2)→(8,3) signal we develop a model based on Einstein rate equations which takes all resonant molecule-electric field interactions into account. We note that at 300 K, states with v = 0 and L = 1-6 are significantly populated, with 27% in L = 2, and 2.4% in states with L ≥ 6. Below, we calculate the Einstein rate coefficients at the level of individual hyperfine states for transitions driven by the laser and BBR fields. Following the approach of Koelemeij [17], we first ignore hyperfine structure and obtain transition dipole moments given by where χ v,L denotes the radial wave function of nuclear motion obtained by numerical solution of the radial Schrödinger equation, R stands for the internuclear separation, and D 1 (R) denotes the permanent electric dipole moment function of the 1sσ electronic ground state.
To take hyperfine structure into account we evaluate the dipole transition matrix element between two states φFS JM J χ v,L and φF S J M J χ v ,L [18,17,19] Equation (5) includes a transformation from the moleculefixed frame to the laboratory-fixed frame [18], with p denoting the polarization state of the electric field with respect to the laboratory-fixed frame (p ∈ {−1, 0, 1}). The linestrength is subsequently calculated by squaring and summing over the M J states [20]: with α ≡ vLFSJ. The Einstein rate coefficients for spontaneous emission, absorption and stimulated emission are subsequently obtained as described in [21], and used in Sec. 4.1 to calculate the expected (0,2)→(8,3) spectrum.

Trapping and cooling Be + and HD +
To achieve narrow linewidths and small systematic shifts, we choose to perform spectroscopy on small samples of HD + molecules in a radiofrequency (rf) ion trap. We reduce the motional temperature of the HD + ions by storing them together with Be + ions which are Dopplercooled by a continuous-wave (cw) 313 nm laser beam (see [17,22,23] for details). The rf trap is placed inside an ultra-high vacuum chamber with a background pressure of 1 × 10 −10 mbar. The rf trap operates at a frequency of 13.2 MHz, leading to Be + radial trap frequencies of ω r = 2π × 290 kHz. The trap geometry and rf circuitry are designed so as to minimize unwanted rf fields and phase differences between the rf electrodes. The two dc electrodes are segmented into two endcaps and a center electrode (Fig. 1). The dc voltages of the center electrodes, rf electrodes and endcap pairs can be individually adjusted to compensate stray electric fields. Be + and HD + are loaded by electron-impact ionization as done by Blythe et al. [24], and monitored with a photomultiplier tube (PMT) and an electron-multiplied charge-coupleddevice (EMCCD) camera. EMCCD images show ellipsoidal mixed-species Coulomb crystals, with a dark core of molecular hydrogen ions surrounded by several shells of fluorescing Be + ions (Fig. 6). The apparatus and procedures for loading and compensation of stray electric fields are described in more detail in [17].  Fig. 1 (Color online) Schematic view of the trap setup. An ultrahigh vacuum (UHV) chamber houses a linear rf trap in which Be + and HD + ions are loaded by electron-impact ionization. Be + ions are Doppler-cooled by the 313 nm laser, which is directed parallel to the trap axis and B-field direction. The 313 nm Be + fluorescence is imaged onto and detected with a PMT and an EMCCD camera. The 782 nm and 532 nm cw REMPD lasers are overlapped with the 313 nm laser and the ions, but propagate in the opposite direction.

Spectroscopy of HD +
The (0,2)→(8,3) transition in HD + is detected destructively through resonance enhanced multi-photon dissociation (REMPD), see Fig. 2a. A 782 nm cw titanium:sapphire laser is used to excite HD + from its vibrational ground state to the v =8 state, which is efficiently dissociated by the field of a co-propagating 532 nm cw laser beam. Both lasers are directed along the trap axis and counter-propagate the 313 nm laser (see Fig. 1). Since all HD + ions are initially in the vibrational ground state (which is too deeply bound to be dissociated by the 532 nm laser), dissociation only takes place if the 782 nm laser is resonant with the (0,2)→(8,3) transition. We used 90 mW of 532 nm radiation, focused to a beam waist of 140 µm, which is sufficient to dissociate an HD + molecule from this level within a few ms, much faster than the spontaneous decay of the v=8 state (lifetime 12 ms). During each measurement cycle, the HD + ions are exposed to the REMPD lasers for 10 s. During the first seconds, the majority of the HD + in L=2 are dissociated, leading to depletion of the L=2 state (Fig. 2). During the remainder of the REMPD period, BBR repopulates the L=2 state from other rotational levels, which enhances the number of dissociated HD + ions by approximately a factor of two.
The 782 nm laser has a linewidth of ∼0.5 MHz and is frequency-stabilized by locking its frequency to a nearby mode of a self-referenced optical frequency comb laser.  To this end, a 63.5 MHz beat note is created by mixing the light of the 782 nm laser and the frequency comb. The frequency comb itself is locked to a rubidium atomic clock for short-term stability, which is disciplined to the 1-pps signal of a GPS receiver for long-term accuracy and traceability to the SI second with 2×10 −12 relative uncertainty.
During REMPD, about 300 mW of 782 nm light is used with a beam waist of ∼120 µm at the location of the ions. The 313 nm cooling laser is detuned by −80 MHz from the cooling resonance and reduced in power to ∼70 µW, which results into an intensity of two times the saturation intensity of the 313 nm cooling transition, and leads to an ion temperature of about 10 mK. To obtain a measure of the fraction of HD + ions lost due to REMPD, we employ so-called secular excitation [25]. This procedure is based on the indirect heating of the Be + which occurs when the motion of the HD + ions is resonantly excited by an additional radial rf field as it is scanned over the resonance frequency. The heating of Be + leads to a change of the 313 nm fluorescence level, which is connected to the number of HD + ions. A single measurement cycle consists of secular scan (10 s), followed by 10 s of REMPD and an-other 10-s secular scan (see Fig. 3) to infer the number of HD + ions remaining after REMPD. To obtain a spectrum, this procedure is repeated for different frequencies of the 782 nm laser (indicated by the variable ν), which are chosen as follows. The ∼25 strongest hyperfine components of the (0, 2) → (8, 3) transition are located in the range (−110 MHz, 140 MHz) around the hyperfine-less frequency. As we expect Doppler broadening to ∼16 MHz, we divide this range into a set of 140 evenly spaced frequencies (spacing ∼ 1.8 MHz) at which REMPD spectroscopy is performed. To convert possible time-varying systematic effects into random noise, we randomize the ordering of the frequency list. For each frequency point, six to seven REMPD measurements are made, which results into a spectrum consisting of 886 REMPD measurements in total.
Scanning the frequency of the radial rf field over the secular motional resonance of HD + (∼830 kHz) temporarily heats up the Coulomb crystal to a few Kelvin. At a detuning of the 313 nm cooling laser of −300 MHz, and using a few mW of laser power (corresponding to ∼ 80 times the saturation intensity), such a 'secular scan' shows up as a peak in the 313 nm beryllium fluorescence versus rf frequency (Fig. 3). In Appendix B we show that the area under this peak, A, scales with the number of HD + ions, although not linearly. We define a spectroscopic signal as the relative difference between the areas of the initial and final secular scan peaks: Repeating this procedure for all 140 pre-selected frequencies of the 782 nm laser while recording fluorescence traces with both the PMT and EMCCD camera, we obtain a spectrum consisting of 1772 data points (ν j , S j ) (with j = 1 . . . 1772).

Spectral lineshape model
In order to determine the hyperfine-less rovibratiotional (0,2)→(8,3) transition frequency, ν 0 , we need a realistic lineshape model which includes the relevant physics present during REMPD. Here the aim is to obtain a spectral lineshape function in which all effects are parameterized. Parameter values are estimated by independent means where possible, or included as a fit parameter otherwise. Importantly, the fit function will contain a variable ν − ν 0,fit , where ν 0,fit is a fit parameter from which we later deduce the value of ν 0 . Before fitting, the (ν, S) data are corrected for reactions with background gas (primariliy H 2 ). This procedure is described in Sec. 4.5.
We start with building a state vector ρ(t) which contains the population in all 62 hyperfine levels in the rotational states ranging from L=0 to L=5 with v=0. This includes 97.6% of the total internal states of v=0 given a BBR temperature of 300 K. We neglect the Zeeman splitting due to 0.19 mT B-field at the location of the trapped ions, as this splitting is negligibly small compared to the Doppler linewidth and the width of the BBR spectrum. The lineshift to ν 0 due to the Zeeman effect is considered separately in Sec. 4.6.2. Also the Stark effect, the electric-quadrupole shift and secondorder Doppler shift are not included in this model, but addressed in Sec. 4.6.2.
During the 10 s of REMPD the hyperfine levels in the L=2 initial state interact with the 782 nm laser and BBR. We here make the simplifying assumption that any population in the v = 8 target state will be dissociated by the 532 nm laser. The interaction with the 782 nm laser is therefore modeled as a simple loss process. The evolution of the state vector ρ(t) is obtained by solving the set of coupled rate equations: where M R and M BBR are the matrices describing the interaction with the REMPD lasers and the BBR field.
Generally, the diagonal elements of these matrices contain depopulating (negative) terms, whereas the off-diagonal elements describe the population a certain state gains from another. Since M R describes losses it only contains diagonal entries of the form where B αα denotes the Einstein coefficient for absorption from a lower hyperfine-rovibrational level α to an upper level α (which for M R the states α and α being restricted to those having v = 0, L = 2 and v = 8, L = 3, respectively). The corresponding transition frequency is ω αα , I laser is the intensity of the 782 nm laser, and D z represents a normalized response function averaged over the distribution of z-velocities of the HD + ions. This involves an integration over all Doppler shifts −k z v z observed by the HD + ensemble, where k z is the wavevector of the 782 nm laser and v z the velocity in the z-direction. If the HD + velocity distribution is thermal (Gaussian), D z depends only on the temperature in the z-direction, T HD + . However, the effects of micromotion and chemistry in the Coulomb crystal during REMPD lead to deviations from a thermal distribution. This implies that D z cannot be described by a single Gaussian lineshape. In Secs. 4.3 and 4.4 these processes are explained in detail and the precise shapes of D z are determined.
The matrix M BBR contains both diagonal and offdiagonal elements, taking into account the rate of exchange of population between all involved levels α and α mediated through all possible electric-dipole transitions, which correspond to spontaneous emission, stimulated emission and absorption by BBR, respectively. The relationship with Eq. (6) is given in [21]. W (ω BBR , T BBR ) denotes the BBR energy distribution function at frequency ω BBR , which is given by: Since the typical frequency of the internal degrees of freedom (> 1 THz) differs from that of the external degrees of freedom (< 1 MHz) by many orders of magnitude, any energy transfer mechanism between them must be of extremely high order and consequently negligibly small. Laser cooling of the external degrees of freedom may therefore be expected to have no significant effect on the temperature of the internal degrees of freedom, which are coupled strongly to (and in equilibrium with) the BBR field [26]. We use Mathematica to solve Eq. (8) in order to obtain ρ(t). We subsequently find the relative HD + loss,  Fig. 3 The secular scan -REMPD -secular scan detection scheme. At the start of each cycle a new sample of HD + ions is loaded into the ion trap. During loading the Be + ions are exposed to neutral HD gas, leading to the formation of a small number of BeH + and BeD + ions (see also Table 1). These are expunged from the trap by applying a dc quadrupole potential of ∼ 0.9 V, which reduces the trap depth such that only ions with mass ≤ 9 amu remain trapped. Five seconds after loading, the HD + is motionally excited by scanning an rf electric field over the secular motional resonance frequency at 830 kHz. A 313 nm laser detuning of −300 MHz is used, in which case the secular excitation results in a rise in the Be + fluorescence. Ten seconds later, the 313 nm laser is tuned to −80 MHz from the Be + cooling transition, and its intensity reduced to I ∼ 2Isat. After 10 s of REMPD, the 313 nm laser settings are restored to their values as used for the first secular scan, and a second secular scan is executed. A smaller 313 nm fluorescence peak indicates loss of HD + .
, by summing over the hyperfine state populations (62 hyperfine states in (v, L) = (0, 2)) before (i.e. t = 0) and after (t = 10) REMPD and computing: Here N i and N f are the numbers of HD + ions present in the trap directly before and after REMPD, respectively. For a thermal ensemble of HD + ions, is a function of the variables ν − ν 0,fit , T HD + , and I laser . In what follows we furthermore assume that T BBR = 300 K, which is the average temperature in our experimental setup.
The question arises what the relation is between the signal S defined in Eq. (7) and the fractional loss defined above. In previous work it was assumed that S and are interchangeable [25,9,27,17]. In Appendix B we study the dependence of the signal S on the initial number of ions N i and the dissociated fraction in detail using realistic molecular dynamics (MD) simulations. We find that the fraction (which is a theoretical construct) can be mapped to the signal S by use of a slightly nonlinear function, whereT 0 is defined as the effective Be + temperature along the z-axis during the initial secular scan of the REMPD cycle (see Appendix B). This means we have to use a five-dimensional fit function An analytical solution of the fit function proves difficult to find, whereas a numerical implementation of the fit function takes excessively long to compute. Therefore we compute values of S fit on a sufficiently dense grid of values ν − ν 0,fit (encompassing 270 values between −120 and 150 MHz), T HD + (20 values between 1 and 20 mK), I laser , andT 0 (9 values between 0.65 × 10 7 and 2.1×10 7 W m −2 ), which we interpolate (again using Mathematica) to find a fast, continuous and smooth approximation to the function , , which we insert into the function S fit . The result is a five-dimensional (5D) fit function which is suitable for nonlinear least squares fitting. A 3D projection of S fit assuming fixed values of I laser ,T 0 and ν 0,fit is plotted in Fig. 4. The reason for treating I laser as a fit parameter instead of inserting a single fixed value is the following. The entire spectroscopy measurement is divided into 15 sessions each taken at a different day. To ensure reproducible laser beam intensities from session to session, we used diaphragms to overlap all laser beams with the 313 nm cooling laser, which itself is aligned with the Be + Coulomb crystal visually using the EMCCD camera. Using a mock version of this setup, we verified that this procedure leads to beam pointing errors up to 40 µm. Assuming a Gaussian distribution of beam pointing errors, we find an intensity at the location of the HD + ions which varies from the intensity in the center of the beam by a factor of 0.6 to 1. Since the spectral line shape is strongly saturated, these variations of the intensity only lead to small signal changes. It is therefore allowed to treat I laser as a free fit parameter which represents the average 782 nm laser intensity for all data points. Similarly, the variables T HD + andT 0 cannot be determined accurately a priori and are treated as free fit parameters as well.

Estimation of absolute trapped ion numbers
As explained in Sec. 4.4, effects of chemistry in the Coulomb crystal significantly influence the measured lineshape of the (0,2)→(8,3) transition. The impact of such effects depends on (and can be estimated from) the absolute numbers of beryllium ions and molecular hydrogen ions. In order to estimate absolute numbers we combine results from MD simulations and spectroscopy. Similar as observed by Blythe et al. [24], our mixed-species ion ensembles contain not only Be + and HD + , but also BeH + , BeD + , H 2 D + , and HD + 2 . In this paragraph we focus on the latter two species which are created during the HD + loading procedure through the exothermic reactions HD + HD + −→ H 2 D + + D (14) and HD + HD After loading, the excess HD gas is pumped out of the vacuum chamber within a few seconds, after which the background vapor resumes its steady-state composition (predominantly H 2 ). Reactions with HD therefore only play a significant role during and just after loading. Triatomic hydrogen ions can be detected by secular excitation. An example is shown in Fig. 5, where the peak in the PMT signal at the left is attributed to the overlapping peaks belonging to species with chargeto-mass ratios 1:4 and 1:5, corresponding to H 2 D + and HD + 2 , respectively. In rf traps, lighter species experience stronger confinement. This is evident from EMCCD camera images which exhibit fluorescing Be + ions surrounding a dark core of lighter species (Fig. 6). The size of the core reflects the total number of light ions. We analyze this by means of MD simulations (Appendix A),  from which a relation is obtained between the size of the dark core and the number of trapped ions. By comparing simulated and real EMCCD images, an average number of ∼ 750 Be + ions is obtained. From the analysis of the image intensity profiles, we cannot distinguish the HD + from the triatomic molecular species. To solve this we use a collection of 140 EMCCD images taken before and after REMPD while the 782 nm laser was tuned at the same fixed frequency near the maximum of the (0,2)→(8,3) spectrum. From the lineshape model we estimate that the average relative HD + loss is 0.57 at this frequency. By comparing the EMCCD images taken before and after REMPD with the simulated images, we also determine the total initial and final numbers of light ions. Combining this with the expected loss of 57 % of the HD + ions, we infer the ratio of HD + numbers to heavier molecular species (H 2 D + , HD + 2 ). An average number of 43 HD + ions is obtained together with 60 ions of heavier species. We cannot determine the relative abundance of H 2 D + and HD + 2 , but previous observations indicate a branching ratio between Eqs. (14) and (15) of 1:1 and, thus, equal abundances [28,29]. The set of EMCCD images shows an appreciable spread in the size of the dark core and, in particular, the ratio of HD + to heavier triatomic hydrogen ions. Variations in both are due to uncontrolled shot-to-shot fluctuations of the HD background pressure during HD + loading. We find somewhat asymmetrical distributions of HD + ions and ions of heavier species, with a wider tail towards higher numbers. The resulting standard deviation of the number of HD + is 41 ions. This also indicates that analysis of EMCDD images (under the present conditions) is not suitable to replace the signal obtained by mass-selective secular excitation (Eq. (7)).
For the treatment of effects of chemistry on the lineshape in Sec. 4.4 below, we consider two scenarios: where scenario b reflects the one sigma upper variation (which well represents the width of the upper tail of the ion number distribution).

Effect of micromotion
The rf quadrupole field of the ion trap inevitably leads to micromotion of ions with non-zero displacement from the trap z-axis. In addition, excess micromotion may be caused by unwanted rf fields arising from geometrical imperfections of the trap electrode structure or phase differences between the rf electrodes [30]. In an ideal linear rf trap, micromotion is strictly radially oriented, but small imperfections in the trap geometry can cause excess micromotion with a component along the trap axis and laser direction, thus adding phase-modulation sidebands to each hyperfine component in the (0,2)→(8,3) spectrum. Due to the combination of an asymmetric and saturated lineshape of the (0,2)→(8,3) spectrum, these sidebands can lead to a shift of ν 0 . Therefore the micromotion amplitude along the 782 nm laser needs to be determined. As the laser propagates virtually parallel to the trap axis, and since the HD + ions are always located near the trap axis, we are primarily concerned with the possible axial micromotion component.
The HD + axial micromotion amplitude x HD + , can be determined through fluorescence measurements of a trapped string of beryllium ions using a modified version of the photon-rf field correlation technique [30]. The idea here is to radially displace a string of about 10 Be + ions by ∼100 µm by applying a static offset field. This will induce significant radial micromotion, in addition to the axial micromotion. The 313 nm cooling laser propagates at an unknown but small angle θ (θ <10 mrad) with respect to the trap axis, and may therefore have a nonnegligible projection along the radial direction. In Appendix C we show that if the rf voltage amplitude, V 0 , is varied, the axial micromotion amplitude scales linearly with V 0 , while the radial one varies as θ/V 2 0 . The latter behavior stems from the V 0 -dependent confinement and the concomitant variation of the Be + radial displacement with V 0 . Thus, measuring the micromotion amplitude for various values of V 0 allows separating the radial and axial contributions.
To determine the micromotion amplitude, we use a similar setup as described in [30]. Photons detected with the PMT are converted to electrical pulses and amplified by an amplifier-discriminator, which generates a START pulse at time t i . Subsequently, a STOP pulse is generated at time t f at the first downward zero crossing of the rf signal. A time-to-amplitude converter (TAC) converts the duration between the START and STOP pulses to a voltage. We record the TAC output voltage with a digital phosphor oscilloscope for 400 ms. We subsequently process the stored voltage trace with a computer algorithm to obtain a histogram of START-STOP time delays, employing 1-ns bins in a range 0-76 ns (i.e. one rf cycle). The bin heights thus reflect the scattering rate as a function of the rf phase, and micromotion will lead to a modulation of the scattering rate about its mean value (Fig. 7). The Be + scattering rate (indicated as R MB , where MB stands for Maxwell-Boltzmann) is given by: where m Be is the Be + mass, T the Be + temperature, I the 313 nm laser intensity, I sat the saturation intensity of the 313 nm cooling transition in Be + , ∆ = 2π×−25 MHz the detuning of the 313 nm laser light, k the wavevector of the 313 nm laser, and v and v µ the secular and micromotion velocities of the Be + ions. While the rf voltage is being varied, the 313 nm laser is displaced so that the ions always are at the maximum of the Gaussian laser intensity profile. The k · v µ term can be written as where x 0,k is the amplitude of Be + micromotion along the direction of the laser wavevector and t 0 is a time offset. We extract x 0,k by fitting Eq. (16) to the acquired fluorescence histogram. Repeating this procedure for various values of V 0 , a list of data points of the form (V 0 , x 0,k ) is obtained. We subsequently extract the radial and axial micromotion components by fitting a model function to these data. The model function is derived in Appendix C.
The procedure of displacing a string of Be + ions and varying V 0 is carried out for both the horizontal and vertical directions. The data and fit functions are shown in Fig. 8 and an average axial micromotion amplitude x HD + (the projection along the 782 nm wavevector) of 11(4) nm is found. As explained in Appendix C, the radial micromotion contribution (due to a possible small angle of the 782 nm laser with the trap axis) averages to zero. We incorporate the micromotion effect by extending the lineshape function D z in Eq. (9) with sidebands at frequencies mΩ with amplitudes J 2 m (k 782 x HD + ). Here, k 782 denotes the wavevector of the 782 nm laser, and the J m are Bessel functions of the first kind, with m an integer in the range [−3, 3].

Effects of chemistry in the Coulomb crystal
During REMPD, H 2 molecules in the background gas can react with the ions in the Coulomb crystal. Such reactions can be divided into two classes: (i) elastic col-lisions, and (ii) inelastic collisions, during which a chemical reaction or charge exchange occurs, and chemical energy is converted into kinetic energy. In general, the kinetic-energy transfer to the ion from elastic collisions with room-temperature particles is much lower than the chemical energy released from inelastic collisions.
At close range r, the electric field of the ion polarizes the neutral molecule which results in an attractive interaction potential U (r) = −αQ 2 /(8π 0 r) 4 . Here α denotes the polarizability volume (in m 3 ) a of the molecule, 0 is the electric constant, and Q is the elementary charge. If we integrate the interaction force over the trajectory of the neutral near the ion (assuming a relatively large impact parameter b > b crit ; see below), we obtain a change of momentum which corresponds to a velocity kick of tens of meters per second for most species.
If a neutral atom or molecule and an ion approach each other within a critical impact parameter b crit = (αq 2 /π 0 µ red v) 1/4 , where µ red and v are the reduced mass and relative velocity of the pair, a so-called Langevin collision occurs, during which the collisional partners spiral towards each other and a chemical reaction can occur at very short range [31]. The chemical reaction products contain hundreds of meV of kinetic energy, which is dissipated into the ion crystal which itself only contains about 2 meV of kinetic energy (at 10 mK). A possible adverse side effect is that such collisions may lead to time-averaged velocity distributions which deviate from thermal distributions. Table 1 shows the relevant reactions during REMPD along with the released chemical energy and their reaction rates. The reaction numbered 1 corresponds to the REMPD process itself. From observations reported in [32] we infer that the ratio of HD + −→ D + + H and HD + −→ H + + D is approximately 1:1. The charge-to-mass ratio of H + is too large for this product to be stably trapped, but the D + ions can stay trapped and can orbit the Coulomb crystal for many seconds. The reaction rate of 1 is calculated from the REMPD model described in Sec. 4.1, and is dependent on the frequency of the 782 nm laser.
Reaction 7 occurs most frequently due to the large number of Be + ions present in the trap. The reaction rate is obtained from the exponential decay of the measured 313 nm fluorescence emitted by a Coulomb crystal of Be + ions, and is in good agreement with the rate estimated from the Langevin cross section given the background pressure of 1 × 10 −8 Pa in our apparatus [33]. The different rate constants of reactions 2 and 3 illustrate the fact that HD + can react with H 2 in two ways: either the H 2 breaks apart, donating an H atom to the HD + molecule, or the HD + breaks apart, after which an H + or D + is added to the neutral molecule. According to [28] and [29] the probability of each scenario is approximately 50 %. In case the ion breaks apart, the a The polarizability volume and the polarizability in SI units are related through α = αSI/(4π 0) Be + ( 2 P 3/2 ) + H2 −→ BeH + + H 0.25 0.0019/0.005 b a The rate (in s −1 per molecule) of D + production is dependent on REMPD time and frequency of the 782 nm laser. b The rate is dependent on the fraction of time a Be + ion spends in the excited 2 P 3/2 state, which is dependent on the 313 nm laser intensity and detuning (−80 MHz or −300 MHz, respectively). v z (m/s) Distribution Fig. 9 The z-velocity distribution of 43 HD + molecules obtained from 100 ms simulation. 14 D + ions with an initial velocity of 6300 ms −1 heat up the ion crystal and give rise to a velocity distribution which differs slightly from a Gaussian (red curve, fitted temperature 11.60(3) mK)) and is better described by a q-Gaussian (green curve, fitted temperature 11.00(3) mK, and fitted q parameter 1.070 (3)).
probability that either a H + or a D + is donated to the H 2 molecule is also 50 %. This leads to a ratio of reaction 2 to 3 of 3:1. The rate of reaction 2 can be measured (keeping in mind that HD + and H + 3 in reaction 3 have the same charge-to-mass ratio and therefore cannot be distinguished) by applying the measurement scheme depicted in Fig. 3 without 782 nm laser, which is further described in Sec. 4.5. The rates of reactions 4, 5 and 6 are obtained from [34]. The kinetic energies of the chemical products are calculated by using the binding energies and energy and momentum conservation laws.
MD simulations show that the fast ionic chemical products may heat up the Coulomb crystal by 1-2 mK (depending on the REMPD rate), and that the HD + velocity distribution becomes slightly non-thermal. A MD simulation of a Coulomb crystal containing 750 lasercooled Be + ions, 40 sympathetically cooled HD + ions and 14 fast D + ions produces the HD + z-velocity distribution shown in Fig. 9. Details of the MD code are given in Appendix A. It turns out that the z-velocity distribution deviates from a Gaussian curve and is better described by a q-Gaussian [35], which is a Gaussian curve with higher wings parameterized by a parameter q. For 1 < q < 3, the q-Gaussian is written as .
(18) Here, ω and ω 0 are the frequency and center frequency, Γ is the gamma function, and β is analogous to the standard deviation of a Gaussian distribution, which is related to the Doppler width. For q = 1, the function reduces to a regular Gaussian distribution. While we use q-Gaussians to describe the simulated data, the observation that q-Gaussians describe the simulated velocity distributions well is merely empirical, and we did not derive this velocity distribution from a particular physical model. Since the value ν 0,fit turns out to be sensitive to the shape of the velocity distribution, it is important that we insert the lineshape based on the correct velocity distribution in Eq. (9), and specify the bounds to within this distribution is valid. An initial analysis reveals that ν 0,fit may shift by several hundreds of kHz by implementing a q-Gaussian with q varying between 1.0 and 1.1. Note that another recent study of MD simulations independently confirmed the non-thermal character of velocity distributions of laser-cooled ion crystals due to collisions with background gas molecules [36].
In the remainder of this section, we determine the q-values applicable to our REMPD measurement with the help of MD simulations. The value of q is dependent on the number of trapped fast ions, N fast . The more fast ions, the higher the q-value. We note that N fast is frequency dependent (more D + are produced near the peak of the REMPD spectrum) as well as time dependent (the production rate of D + is governed by the rate equations, Eq. (8)).
When a fast ion collides with a cold ion, each particle may undergo a nonadiabatic transition to a different solution of the Mathieu equation which governs its motion in the trap [37]. This implies that a fast ion can either lose or gain energy during such a collision. Since the energy change per collision is relatively small, fast ions can retain their high speeds in the trap for many seconds. Values for N fast can be obtained by solving the rate equation: The term α prod is the production rate from a number of N source particles, such as Be + or HD + . In Sec. 4.2 we determined N source,Be + =750 and N source,HD + = 43 or N source,HD + = 84. α relax is the rate at which fast ions are cooled and become embedded within the Coulomb crystal, while α esc is the rate at which (fast) ions escape from the trap. The values of α prod correspond to the reaction rates in Table 1. However, obtaining α relax and α esc requires a multitude of individual MD simulations, with simulation periods of several seconds each. Even with current academic supercomputers, the total time to perform such simulations is prohibitively long. Therefore, we consider two extreme scenarios for relaxation and escape rates of trapped fast ions: -Scenario 1: A minimum number of fast ions is present in the trap. α relax and α esc are set to their maximum value of one per second for all species, which is based on the observation that no ion loss and no relaxation are observed over simulated times up to 800 ms. This scenario will produce the smallest value of q. -Scenario 2: A maximum number of fast ions is present in the trap, which is realized by setting α relax and α esc to their minimum value of zero. All fast ions remain in the trap at high speed for the entire 10 seconds of REMPD. This scenario leads to the largest value of q.
We make the assumption that fast ions do not mutually interact when present in numbers of ten or less, so that the observations based on MD simulations with ten fast ions are also valid for smaller numbers of fast ions. Note that BeH + ions are already created during the first secular scan before the REMPD phase starts. Fast specimen of H 2 D + and H + 3 occur less frequently in the trap, and in line with the extreme scenarios introduced above we assume zero H 2 D + and H + 3 for scenario 1, and 3 H 2 D + and 1 H + 3 for scenario 2. Together with scenarios a and b described in Sec. 4.2 this gives us four scenarios in total, and therefore four different spectral fit functions and four different ν 0,fit results.
For all possible combinations of fast ions present during REMPD (e.g. 1 D + , 2 BeH + and 1 H 2 D + or 3 D + , 3 BeH + and 3 H 2 D + ) an MD simulation is carried out. From these simulations, the HD + z-velocity distribution is determined and a q-Gaussian is fitted which results in one q value for each simulated case. As mentioned before, the production rate of fast D + depends on the REMPD rate (and thus on the 782 nm laser frequency ν), which itself depends on the time-dependent number of available HD + ions in the target state. To take this properly into account we introduce a time and frequency dependent parameter q(t, ν) as follows. For each of the four scenarios, the number of fast D + is simulated on a grid of different REMPD durations, t j , and of different frequencies, ν j of the 782 nm laser. On each point of this two-dimensional grid, the number of fast D + is combined with the number of other fast ions, and the corresponding value of q(t j , ν j ) is looked up in a library of q values, obtained from many MD simulations performed with various combinations and abundances of fast ion species. Interpolation of the grid q(t j , ν j ) leads to a smooth continuous function q(t, ν), which is subsequently inserted into the lineshape function D z used in Eqs. (8) and (9). Figure 10 shows the 3D plots of q(t, ν) for the different scenarios.
Besides the q value, also the ion temperature T HD + is frequency and time dependent. Due to a larger number of D + ions at the top of the spectrum than at the wings, the temperature differences between top and wings can reach a few mK. We note that the increase of q and T HD + share the same origin (namely collisions with fast ions), and in Fig. 11 we show the relation between q and T HD + , obtained from fitting q-Gaussians to simulated velocity distributions. We find a linear dependence of the form T HD + [1 + h(q(t, ν) − 1)], with the slope h obtained from the fit. In scenario 1, the temperature difference between top and wings is found to be 0.5 mK. For scenario 2, the q Fig. 11 The relation between the temperature of a simulated Coulomb crystal and the q value. Each point represents a simulation of approximately 100 ms. The data are best fitted with a linear function (blue line). estimated temperature difference is 2.5(5) mK, where the uncertainty of 0.5 mK is treated as one standard deviation. The t and ν dependent temperature is also included in D z .
A schematic overview of the various steps involved in the construction of the fit function S fit , as well as the role of the four scenarios, is presented in Fig. 12.

Background gas reactions
During the REMPD phase the number of HD + ions is not only reduced through photodissociation by the lasers. As described in Sec. 4.4, trapped ions can react with residual H 2 molecules in the vacuum setup. In order to correct the spectroscopic signal for these so-called background losses, the rates of reactions 2 and 3 in Table 1 are measured and included into the rate equations, Eq. (8), as an additional loss channel.
The spectral data were acquired over the course of several months during 15 independent measurement sessions lasting several hours each (Fig. 13). During this period the background pressure varied from session to session. For each session, the HD + background signal is obtained by using the measurement scheme depicted in Fig. 3, but with a shutter blocking the 782 nm laser, thus preventing REMPD. Typically, a few HD + ions react with H 2 which is detected as a small difference between the secular scan peak areas A bg,i and A bg,f . During each measurement session, a series of ∼ 7 background loss measurements is carried out three times, providing a data set of about 20 background measurement per session. From the average background loss signal per session a reaction rate γ bg is extracted from the relation Here, the vector ρ bg equals ρ of Eq. (8) extended by two additional rows which describe the occurence of ions in the form of H 2 D + or H + 3 . M bg is a diagonal matrix containing γ bg which describes the HD + losses. We take into account the fact that conversion of HD + to H + 3 in reaction 3 in Table 1 is not detected by the method of detection through secular excitation because HD + and H + 3 have the same mass-to-charge ratio. This also means that the measured values of γ bg only represent the rates for reaction 2. As explained in Sec. 4.4, the reaction rate of 3 is a factor of three lower, which is taken into account as well.
We correct the data set of each session individually for the background signal following an iterative procedure. During the first step we insert an initial (coarse) estimate of the value ofT 0 in Eq. (20), and we simply subtract the signal predicted by the model without background losses (based on Eq. (8)) from the signal prediction including background losses (based on Eq. (21)). In this way an estimate of the background signal is obtained which is subsequently subtracted from the raw measurement data. Then the spectral fit function S fit (see Sec. 4.1) is fitted to the corrected data points with free fit parametersT 0 , T HD + , I laser and ν 0,fit . To find an improved estimate of the background signal, the thus found values ofT 0 , T HD + , I laser are reinserted into Eqs. (8) and (21), and the above procedure is repeated. After a few iterations, convergence is achieved.
The apparent loss of HD + ions is also influenced by reactions 4 and 5 in Table 1, which result in the production of additional H + 3 during REMPD, leading to a reduction of the observed REMPD and background signals mentioned above. This effect cannot be easily assessed experimentally, and instead we estimate it for both scenarios a and b using the reaction rates of Table 1, leading to slightly modified values of N i and N f . This translates to a small correction to in Eqs. (11) and (20); see also Fig. 12). In addition, reactions with background gas also take place during the secular scans, thus influencing the determination of the initial and final numbers of particles with charge-to-mass ratio 1:3 itself. For example, the conversion of HD + to H 2 D + (reaction 2) and to H + 3 (reaction 3) have a small effect on the number of estimated HD + ions from the secular scan. Again, we can readily correct for these effects, knowing the values of all the relevant reaction rates, and noting that the reactions take place on a much longer timescale than the secular scan itself so that constant production rates can be assumed. Altogether, the corrections to N i and N f are at the level of a few percent, and are applied through Eqs. (11) and (20) (see also Fig. 12).

Spectrum, systematic effects and final result
The function S fit (Eq. (13)) is fitted to the REMPD data set after correction for the background signal. Figure 14 shows the REMPD data set and fit function for scenario 1a. The noise in the spectrum has several origins. Firstly, the number of trapped ions is relatively small and varies from shot to shot. Secondly, the population in the various hyperfine states of the L=2 state varies from shot to shot, as expected for hyperfine states with a mean occupancy of order unity. The (stochastic) BBR interaction, which couples the states with L = 2 with other rotational states, introduces additional random signal variations. Furthermore, the variation of the number of reactions of HD + with the background gas is in the order of a few per shot, which dominates the noise for low REMPD signals. Finally, part of the noise originates from random intensity variations due to spatial alignment variations of the 313 nm, 782 nm and 532 nm lasers.
For each of the four scenarios (1a,1b,2a,2b) we obtain a particular set of fit parameters, which are listed in Table 2. The correlation coefficients of the fit parameters are presented in Table 3 and the values for ν 0,fit −ν th are graphically shown in Fig. 15. The error bars represent the ±1σ fit uncertainty which can be considered as the purely statistical precision of the spectroscopy measurement. We remark that the sensitivity to the chemistry processes (scenarios 1 and 2) is much stronger than the sensitivity to the numbers of HD + ions in the trap (scenarios a and b), and that the values ν 0,fit for 1a and 2a represent extreme upper and lower limits (with respect to the line shift due to chemistry) for the average number of HD + ions. We therefore chose to obtain our final result for ν 0,fit by taking the mean of these two values, while in- terpreting the mean single-fit error of (which is virtually the same for all four results) as the statistical uncertainty of the final result. We subsequently quantify the 'whichscenario' uncertainty as follows. For scenarios 1a and 1b (and similarly for 2a and 2b) the difference is due to a 1σ variation in the number of HD + ions. Therefore, we treat the frequency interval between the two values corresponding to scenario a and b as the corresponding 1σ interval, which amounts to 80 kHz when averaged over the two scenario's 1 and 2. To find the error corresponding to scenarios 1 and 2, we take the frequency interval between the values found for scenarios 1a and 2a (which are essentially extreme limits) and conservatively equate the interval to a 68 % confidence interval. The interval thus corresponds to 2σ, with σ = 0.23 MHz. The uncertainty of 0.5 mK in the temperature difference between top and wings of the spectrum (see Sec. 4.4) results in 28 kHz difference in ν 0,fit , which is treated as a 1σ variation. The frequency shifts due to these systematic effects are listed together with their uncertainties in Table 4.

Frequency uncertainty of the 782 nm laser
The beat note of the frequency-locked 782 nm laser with the optical frequency comb is counted during REMPD. We use the beat-note frequencies to compute the Allan deviation, which is of the order of 0.1 MHz after 10 s averaging. The uncertainty of the 782 nm laser frequency can be transferred to an uncertainty in the (0,2)→(8,3) fit result ν 0,fit by taking the Allan deviation as a measure of the standard deviation of a Gaussian noise distribution, describing the laser frequency offset from the set frequency during each REMPD cycle. We perform a Monte Carlo simulation in which each of the 140 measurement frequencies is assigned a frequency offset, selected at random from the Gaussian distribution. Repeating this 100 times generates 100 different spectral data sets. Fitting S fit to each of the data sets gives 100 different values of ν 0,fit . A histogram of the resulting distribution of ν 0,fit values is shown in Fig. 16. From the histogram we find a mean offset 0.5 kHz from the frequency value ν 0,fit found for scenario 1a, and a standard deviation of 8 kHz. We conclude that frequency noise introduces no significant bias, and we conservatively assume the uncertainty of the frequency measurement to be 0.01 MHz. The 782 nm laser has a Gaussian lineshape with a width of ∼0.5 MHz. The convolution of this lineshape with the (Gaussian) Doppler-broadened line (16 MHz) will give rise to another Gaussian lineshape. Since the linewidths add up quadratically, the increase in linewidth is smaller then the uncertainty of the linewidth due to the fit uncertainty of the temperatures, which is ∼ 0.8 mK. Therefore, we consider the laser linewidth to be completely absorbed into the fitted temperature T HD + with no significant effect on its value. 4.6.2 Zeeman, Stark and other shifts So far we have neglected the Zeeman splitting of the lines in the spectrum. Incorporating the Zeeman effect makes the hyperfine transition matrices very large and Mathematica is only able to solve the rate equations [Eq. (8)] effectively if lineshapes are not too complicated. We circumvent these issues as follows. First, we calculate the lineshifts and linestrenghts of the magnetic subcompo- Table 3 The correlation coefficients for the fit of S fit (scenario 1a) to the data. This procedure is similar to that followed in Refs. [9,38], and can readily be done for the case of σ + , σ − and π transitions, and superpositions thereof. As the Zeeman splitting is small compared to the Doppler width, the magnetic subcomponents belonging to the same hyperfine component overlap well within the profile of the lineshape function D z , forming a new composite (and Zeeman-shifted) lineshape function D z . This new lineshape function is subsequently used in Eq. (9). For simplicity, we do not implement effects of micromotion and chemistry in this analysis, and we compare ν 0,fit fit results based on versions with B = 0.19 mT with a version with zero B field. It is important to note that the linear polarization of the 782 nm laser is practically perpendicular to the B field, so that during the (0,2)→(8,3) excitation only σ + and σ − transitions are driven. Due to polarization imperfections, caused by the polarization optics and the slightly birefringent viewports of the vacuum chamber, the two circularly polarized components are estimated to have a maximum possible intensity imbalance of 2%. Figure 17 shows the offset between the ν 0,fit fit results, obtained with the above model for several σ + /σ − intensity ratios, and the fit result we found previously assuming zero magnetic field. gravity of the spectral line and can therefore be considered as weighted means of all the shifts of the single hyperfine components. The calculation of the Stark shifts is shortly explained in Appendix D. Stark shifts due to the BBR and trap rf field are calculated in [19] and are smaller than 1 Hz. Together, this gives us a total Stark shift of 1.3(1) kHz. The uncertainty in this value stems almost exclusively from the accuracy to which the laser beam intensities are known.
A conservative upper limit of 100 Hz to the electricquadrupole shift for the (v,L):(0,2)→(8,3) transition is obtained from [39]. The second-order Doppler effect is dominated by average micromotion velocity of the HD + ions, and is found to be less than 5 Hz. The values of the aforementioned shifts and their uncertainties are listed in Table 4.
If we compare ν 0,fit from scenario 1a based on the model containing rotational states L ≤5 with an extended version containing L ≤6 states we find a shift of ν 0,fit of 28 kHz, which we treat as the uncertainty due to the neglect of population in L = 6. The rateequation model furthermore includes the BBR temperature, which we estimate to be 300 K with an uncertainty of about 5 K, caused by day-to-day variations of the temperature in the laboratory, and by a possibly elevated temperature of the trap electrodes due to rf current dissipation. If we compare the ν 0,fit values after inserting BBR temperatures of 300 K and 305 K, we obtain a difference of 5 kHz, which we include in the uncertainty budget. Micromotion also causes a shift: by comparing the ν 0,fit value from scenario 1a that includes micromotion (amplitude 11(4) nm) with a the result of a version with zero micromotion, a shift of 0.055 (20) MHz is obtained. Finally, we investigated the effect due to the uncertainty of the spin coefficient E 4 (see Eq. (1)), which is estimated to be 50 kHz [16]. Comparing fits with E 4 values differing by 50 kHz we find that this has a negligibly small effect of ∼1 Hz on the result for ν 0 . 1 a This is a shift with respect to a scenario with a zero effect of q or micromotion, and serves to illustrate the size of the effect. The shift itself, however, is absorbed in the value of ν 0,fit . b The value of these shifts is actually nonzero but negligibly small, and therefore ignored here.  Table 4 shows the error budget, with a total uncertainty of 0.41 MHz that corresponds to 1.1 ppb. This result differs 0.21 MHz (0.6 ppb) from the more accurate theoretical value, ν th = 383,407,177.150(15) MHz [2]. The two main uncertainty contributions are the statistical fit error of 0.33 MHz and the uncertainty in the q-factor scenario of 0.23 MHz.

Conclusion, implications and outlook
We have measured the (v,L):(0,2)→(8,3) transition in the HD + molecule with 0.85 p.p.b (0.33 MHz) precision, which is the first sub-p.p.b. resolution achieved in molecular spectroscopy. A thorough analysis of systematic effects points out that the total uncertainty is 1.1 p.p.b., and the result (ν 0 = 383, 407, 177.38(41) MHz) differs by only 0.6 p.p.b. from the theoretically predicted value (ν th = 383, 407, 177.150(15) MHz). A large contribution to the systematic uncertainty is the effect of chemical reactions in the Coulomb crystal, of which the 1σ uncertainty is 0.61 p.p.b (0.23 MHz). This effect, which had not been recognized before, causes a nonthermal velocity distribution that can be approximated by a q-Gaussian function, and which significantly influences the accuracy of laser spectroscopy of composite lineshapes in the presence of strong saturation and depletion of the HD + sample. This is a situation regularly encountered in laser spectroscopy of Doppler broadened transitions in finite samples of trapped molecular ions.
The agreement between experimental and theoretical data has several implications. First, it enables a test of molecular theory at the 1-p.p.b. level, and a test of molecular QED at the level of 2.7 × 10 −4 , which are the most stringent tests performed so far. Second, it allows us to put new bounds on the existence of hypothetical fifth forces, and put new limits on the compactification radius of higher dimensions, as described in detail in [3,4,11]. Third, the result presented here can be used to obtain a new value the proton-to-electron mass ratio with a precision of 2.9 p.p.b. [11], for the first time from molecular spectroscopy as proposed already four decades ago [5].
Our analysis clearly demonstrates that the first-order Doppler effect is responsible for the largest contribution to the uncertainty. In fact, removing the first-order Doppler effect would render the frequency measurement as the largest source of error (0.026 p.p.b), thereby immediately improving the uncertainty by about a factor of 40. Such -and even larger-improvements are possible using two-photon spectroscopy. For example, in Refs. [21,40] an experiment was proposed in which the (v, L) = (0, 3) → (9, 3) line in HD + is addressed through a twophoton transition with nearly degenerate photons. Using counter-propagating laser beams with a narrow linewidth, the Lamb-Dicke regime may be reached in the present apparatus, such that first-order Doppler broadening is entirely eliminated, while all other systematic effects could be controlled below the 1×10 −13 uncertainty level. Such spectroscopy would allow more stringent tests of molecular theory and QED, tighter bounds on new physics at theÅngström scale, a competitive determination of the proton-electron mass ratio, and even contribute to the determination of several other fundamental constants including the Rydberg constant, the deuteron-electron mass ratio, and the proton and deuteron charge radii [40].

Acknowledgements
We are indebted to J. Bouma, T. Pinkert and R. Kortekaas for technical assistance, and to V. Korobov, E. Hudson and R. Gerritsma for fruitful discussions. This research was funded through the Netherlands Foundation for Fundamental Research on Matter (FOM), the COST action MP1001 IOTA, and the Dutch-French bilateral Van Gogh Programme. J.C.J.K. thanks the Netherlands Organisation for Scientific Research (NWO) and the Netherlands Technology Foundation (STW) for support. SURFsara (www.surfsara.nl) is acknowledged for the support in using the Lisa Compute Cluster for MD simulations.

Appendices
A Molecular Dynamics simulations MD simulations are implemented in Fortran code in order to realistically describe the dynamics of trapped and laser-cooled ions in the presence of the time-dependent trapping field, 313 nm photon scattering by the Be + ions, and fast ionic products from chemical reactions. A cycle of one time step starts by computing the sum of the forces acting on each ion, which consists of the Coulomb force, F C , the time-dependent force from the trapping field field, F trap , and an optional rf electric field, F SS , which drives the secular motion: The radial and axial part of F trap are given by and F trap,z = az + bz 3 + cz 5 , where ω z is the secular angular frequency in the z-direction. The constants a, b and c depend on the trap geometry, which are determined through a finite-element analysis performed with the software package SIMION. The forces exerted on each ion are calculated, and trajectories are obtained using the Velocity Verlet method [41]. Doppler-cooling is included at the level of single-photon scattering. Photon momentum kicks are simulated as velocity changes where absorption only takes place in the laser direction. In order to include ion motional heating which occurs in the trap, we implemented additional stochastic velocity kicks with a size of the recoil momentum of a single 313 nm photon with random directions. If an average kick rate of 75 MHz is used, ion temperatures of around 10 mK are obtained. The processes of elastic and inelastic neutral-ion collisions are simulated as velocity kicks in random directions. For example, simulating reaction 7 in Table 1, a Be + ion is substituted with a BeH + ion at 10 mK, after which its speed is modified so as to give it 0.25 eV of kinetic energy.
Simulation of particles which in some cases have high velocities requires the use of a variable time step size, ∆τ . If the proper step size is not observed, two particles with a high velocity difference at close distance could 'skip' each other within one time step instead of colliding. The default step size is ∆τ = 0.2 ns. However, if for any of the trapped particles the condition v j ∆τ > 10 × min{∆x jk } is met, where min{∆x jk } is the distance between particle j with velocity v j and the nearest particle k, the time step ∆τ is reduced by a factor of 10. Likewise, the step size is increased if the colliding particles separate again and v j ∆τ < 100 × min{∆x jk }.
To simulate EMCCD images, we made use of a simpler MD implementation, which treats the motion of the ions in the pseudopotential approximation and which does not include high-energy ions. This allows for a larger integration time step (10 ns) and, thus, faster MD simulations.

B Derivation of the non-linear fluorescence function
The relative HD + loss during REMPD, , is related to the spectroscopic signal S through the non-linear function f NL , which we derive here. The fluorescence yield during a secular scan depends on the Be + temperature, T , and is described by the scattering rate formula integrated over a Maxwell-Bolzmann velocity distribution. Neglecting micromotion effects, we take the scattering rate R MB = R MB (T, ∆, I/I sat ) defined in Eq. (16). During a secular scan in the experiment we use the values ∆ = 2π × −300 MHz and I/I sat = 67, and in what follows we drop these variables from the function argument of R M B . While performing the secular scan, the temperature T varies, which leads to a fluorescence peak as described by Eq. (16). The spectroscopic signal S is the relative difference between the areas under the fluorescence peaks (see Eq. (7)). We may rewrite the area, A, as where C is a constant taking into account the collection and quantum efficiencies of the PMT or EMCCD imaging system, ∆t denotes the duration of the secular scan (10 s), andT stands for the 'effective' value of T during a secular scan. The effective temperatureT is defined through Eq. (25), and can be interpreted as follows. If we would fix the secular excitation field frequency and amplitude at a certain value during a scan, the Be + temperature and the fluorescence level would remain constant.T is the constant temperature that leads to the same area under the fluorescence trace as during a true secular scan. Likewise,T bl stands for the effective baseline temperature during ∆t (corresponding to the fluo-  Fig. 18 A simulated (yellow) and a real (blue) secular scan peak plotted on the same frequency axis. The signals are scaled vertically to achieve matching peak heights.
rescence level that results if no HD + ions are present). In the experiment, the baseline may have a small slope due to the wing of the secular resonance of particles with mass 4 and 5 amu (see, for example, Fig. 5). This slope is detected and removed by the Mathematica code we use to analyze the PMT signal traces.
Both in the experiment and in the MD simulations described below, we observe the area A under experimental or simulated fluorescence traces, to which we subsequently can assign an effective temperatureT through Eq. (25). We emphasize that no attempt is made to deriveT directly from, for example, the simulated velocities of Be + ions. Also note that in practice we only use Eq. (25) to assign effective temperatures to simulated fluorescence traces, as in this case the constant C is known (C = 1).
Inserting Eq. (25) into Eq. (7), we obtain which is a relationship between the spectroscopic signal S and the effective ion temperatures during the initial and final secular scan,T i andT f respectively. The relationship between the number of trapped HD + molecules andT can be obtained from MD simulations. A Coulomb crystal is simulated containing 750 trapped Doppler-cooled Be + ions and with HD + numbers varying from 0 to 100. This is done once using a number of additional H 2 D + and HD + 2 ions equal to that of scenario a, and once using the numbers of H 2 D + and HD + 2 of scenario b (see Sec. 4.2). The simulated secular scans over the HD + secular resonance frequency produce fluorescence peaks which agree qualitatively with those obtained in the laboratory as shown in Fig. 18. In the experiment, secular scans are acquired over a time span of 10 s. The dynamics due to the time-varying frequency of the ac electric field take place at a time scale much longer than the time scale of fluorescence dynamics during laser cooling, which takes place at time scales of the order of 10 µs [42], and also longer than typical ion oscillation periods and the frequency of the ac electric field itself (1-20 µs). However, due to limited computational resources, the simulated duration of a secular scan is approximately 100 ms, which implies that the simulated dynamics take place at a considerably faster rate than in the experiment, typically on the scale of milliseconds. However, this still is much longer than the time scale of fluorescence dynamics and the motional dynamics. Therefore, we assume that the simulated secular scan peaks provide a reliable model of the experimentally observed secular scans.
The MD simulations reveal a linear relationship between the number of trapped HD + ions andT . This agrees with the intuitive picture of Be + ions with frictionally damped motion (because of the laser cooling), whose temperature rise during secular excitation is directly proportional to the number of HD + ions. Figure 19 shows the (N HD + ,T ) relationship for the two scenarios a and b. Having established thatT is a linear measure of the number of trapped HD + molecules, we now combine the relations = (N i − N f )/N i ,T i = c 1 N i + c 2 and T f = c 1 N f + c 2 , where c 1 and c 2 are constants derived from MD simulations (Fig. 19), to obtain T i can also be defined as the effective temperature with zero HD + loss (T i =T f ( = 0) ≡T 0 ), while the term c 2 can be considered as the effective baseline temperature (c 2 =T f ( =1)=T bl ). Inserting Eq. (27) into Eq. (26) results in the nonlinear function which is plotted for scenario a in Fig. 20. Note that the nonlinear dependence on originates from the nonlinear dependence of R MB on T in Eq. (16).
In the analysisT 0 is treated as a free fit parameter.T bl is kept at a fixed value which is obtained from MD simulations. From Fig. 19 it can be seen thatT bl 0.5 K for scenarios a and b. This indicates that the rf field used for secular excitation already induces heating of Be + while the field is still far away from the Be + resonance (at ∼ 300 kHz). This effect is also seen in the experiment.
The nonlinear function f NL is used to map the relative HD + loss onto the spectroscopic signal S. However, for the correction of background signals (Sec. 4.5) we need to map S to , which requires the inverse nonlinear function f −1 NL . This inverse function is obtained numerically by use of Mathematica.

C Micromotion fit function
The time-dependent electric field of the trap, E t , can be expressed as [30] E t (x, y, z, t) ∼ = − V 0 R 2 (xx − yŷ) cos(Ωt) where R is half the distance between two diagonally opposing electrodes, U 0 is the endcap voltage, Z 0 stands for half the distance between the end caps, and κ is a shielding factor. The Be + micromotion amplitude can be written as The measured micromotion amplitude can be written as which is the projection of the 313 nm laser direction onto x 0 . From simulations of the rf trap circuitry with the simulation software Spice, we find a small possible phase difference φ ac of 4 mrad in between the rf electrodes, which has a negligible effect on the ion micromotion and is ignored here. Using the program SIMION, we calculate the shielding factor κ and the static electric fields E dc as a function of the dc voltages applied to the trap electrodes. From the static electric fields, the radial ion displacement r d is obtained by balancing the ponderemotive force and static E-field in the radial direction, where ω r is the radial secular trap frequency. By inserting the x-and y-components of r d into Eq. (29), the value of E t at the location of the ions is obtained. A geometric imperfection of the trap could lead to an axial rf field, which can be written as (here we ignore the small modification of the radial rf field of the trap due to the same imperfection): where x HD is the HD + micromotion amplitude along the trap z-axis, m HD is the mass of HD + , and V 0,e is the rf voltage used during the spectroscopic measurements, which is 270 V. Now, we turn to the case of a linear string of Be + ions, which is the configuration used to determine the axial rf field amplitude. Adding E ax,HD + to the z-component of E t gives a new expression for E t which is inserted into Eq. (30). We then obtain the following expression for x 0 : which is subsequently inserted into Eq. (31), together with the wavevector, which is written as k = 2π λ (sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)) .
Here θ is the angle between k and the trap z-axis, and φ is the angle between k and the trap y-axis, which is very close to π/4 in our setup. The value of θ lies between ±10 mrad and is treated as a free fit parameter. We insert Eqs. (35) and (34) into Eq. (31) and then expand the expression in powers of θ. This gives us the following fit function: Here E h,offs , E v,offs are the applied static electric fields in the horizontal and vertical directions, respectively, and δE h , δE v are the unknown offset electric fields (due to e.g. charging of electrodes). The Mathieu parameters a M and q M are given by The displacement of the Be + string in the vertical direction can be accurately determined with images of the EMCCD camera, and therefore δE h can be zeroed (for example by minimizing the displacement of the Be + string while the radial confinement of the trap is modulated by varying the rf amplitude). However, the displacement in the horizontal direction (i.e. perpendicular to the EMCCD image plane) is not accurately known and therefore we treat δE h as another free fit parameter. In summary, we use Eq. (36) as a fit function with x HD + , θ and δE h as free fit parameters while neglecting higher orders of θ. The fitted curves and the result for x 0,k are shown in Sec. 4.3.
The question arises what happens if the 782 nm laser propagates at a small angle with respect to the trap axis, while the HD + ions form a shell structure around the trap axis. In this case a small fraction of the radial micromotion is projected onto the wavevector. However, from Eqs.  it follows that the sign of this additional micromotion alternates for each quadrant in the (x, y) plane. As long as the radial micromotion component does not exceed the axial micromotion amplitude (which is the case here), the former averages out to zero given the radial symmetry of the HD + crystal.

D Stark shift calculations
Here we summarize the formulas that are used to calculate the ac Stark shift of a rovibrational transition (v, L) → (v , L ) in the HD + molecule induced by a laser with intensity I and polarization state p. A general expression for the second-order energy shift depending on the angle θ between the polarization direction and the quantization axis is: where P 2 (x) = 1 2 (3x 2 −1) is a Legendre polynomial. This expression contains the scalar and tensor polarizabilities where a 0 is the Bohr radius and Q s and Q t stand for the two-photon scalar and tensor matrix elements: Here Q (0) and Q (2) are the irreducible scalar and tensor components that belong to the two-photon operator (in atomic units): with Hamiltonian H, dipole moment operator, d and polarization vector p . The matrix elements Q s and Q t were calculated numerically using the three-body variational wave functions described in [8].
Since the hyperfine structure is partially resolved in this spectrum, we also have to consider the contribution of the Stark shifts to off-resonant coupling to hyperfine levels in v = 0 and v = 8 by the 782 nm laser during spectroscopy. Here, the situation is more complicated as the 782 nm laser also non-resonantly couples v = 8 states to continuum states above the dissociation limit of the 1sσ electronic ground state. The 782-nm contribution to the Stark shift was calculated at the hyperfine level, and will be published elsewhere. The corresponding shifts turn out to be negligible for our experiment, contributing only at the level of a few Hertz.