Highprecision spectroscopy of the HD molecule at the 1p.p.b. level
Abstract
Recently we reported a high precision optical frequency measurement of the : (0,2)(8,3) vibrational overtone transition in trapped deuterated molecular hydrogen (HD) ions at 10 mK temperature. Achieving a resolution of 0.85 partsperbillion (p.p.b.) we found the experimental value ( MHz) to be in agreement with the value from molecular theory ( 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 protonelectron 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 (laserinduced) chemical reactions, as such deviations turn out to significantly shift the hyperfineless vibrational frequency as inferred from the saturated and Dopplerbroadened spectrum, which contains partly unresolved hyperfine structure.
1 Introduction
Because of their relative simplicity, molecular hydrogen ions such as H and HD are benchmark systems for molecular theory Korobov2014a (); Korobov2014b (), and suitable probes of fundamental physical models Salumbides2013 (); Salumbides2015 (). Already four decades ago, Wing et al. Wing1976 () performed measurements of rovibrational transitions in HD and conjectured that such measurements could be used to test quantum electrodynamics (QED) theory in molecules, and lead to an improved value of the protontoelectron mass ratio, . Today, QED corrections to rovibrational transition frequencies in H and HD have been calculated up to the order , with the electron mass and the finestructure constant, leading to relative uncertainties of the order of – Korobov2014a (); Korobov2014b (). The frequencies of these transitions are typically in the (near) infrared with linewidths below 10 Hz, which makes them amenable to highresolution laser spectroscopy. Optical clocks based on trapped atomic ions have shown that laser spectroscopy at very high accuracy (below one part in ) is possible Chou2010 (). Recent theoretical studies point out that also for molecular hydrogen ions experimental uncertainties in the range should be possible Schiller2014 (); Karr2014 ().
Recent frequency measurements of rovibrational transitions in HD were done with a reported uncertainty of 1–2 p.p.b. Koelemeij2007a (); Bressel2012 (). 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 Bressel2012 (). 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 highprecision frequency measurement of the rovibrational transition (v,L):(0,2)(8,3) in the HD molecule, and compare it with stateoftheart 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 Biesheuvel2015 ().
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 lasercooled beryllium ions. The analysis of the spectroscopic data, systematic frequency shifts and the results are discussed in Sec. 4, followed by the conclusions in Sec. 5.
2 Theory
2.1 Calculation of rovibrational frequency transitions in HD
The calculation of rovibrational energies in quantum mechanical threebody systems is usually divided into a nonrelativistic part, complemented by the relativistic and radiative parts, and finite nuclear size contributions. The nonrelativistic part is calculated through solving the threebody Schrödinger equation, which can be done with practically infinite precision Korobov2000 () (up to a relative precision of Li2007 (); Hijikata2009 ()). The resulting wavefunctions allow an analytical evaluation of the BreitPauli Hamiltonian and the leadingorder radiative corrections.
Recently, the precision of relativistic and radiative energy corrections to the nonrelativistic energies was strongly improved. With the inclusion of the full set of contributions of order and leadingorder terms of order , the relative uncertainty is now below . For example, the theoretically determined value of the HD (v,L):(0,2)(8,3) transition frequency, , is 383,407,177.150(15) MHz, and has a relative uncertainty of Korobov2014a (); Korobov2014b (). 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 uncertainty of the 2010 Committee on Data for Science and Technology (CODATA) value of CODATA2010 (), which translates to a frequency uncertainty of 59 kHz.
2.2 Hyperfine structure and rotational states
Since the HD constituents possess nonzero spin, the rovibrational transition spectra contain hyperfine structure due to spinspin and spinorbit couplings. In Bakalov2006 () the hyperfine energy levels and eigenstates in HD are calculated by diagonalization of the effective spin Hamiltonian
(1)  
where the spin coefficients, , are obtained by averaging the BreitPauli Hamiltonian over the nonrelativistic wavefunctions, is the rotational angular momentum operator, and , and are the electron, proton and deuteron spin operators. and are spherical tensors composed of angular momenta, whose explicit form is given in Bakalov2006 (). The strongest coupling is the protonelectron spinspin interaction (the term in in Eq. (1)), and the preferred angular momentum coupling scheme is
(2) 
For states with , and 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 , with the magnetic quantum number corresponding to the projection of onto the laboratoryfixed zaxis. Note that after diagonalization the quantum numbers are only approximately good, and a hyperfine eigenstate can be expressed in the ’pure’ basis states in the form
(3) 
with realvalued coefficients . In Bakalov2006 () 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 hyperfineless (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).
2.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 moleculeelectric field interactions into account. We note that at 300 K, states with and –6 are significantly populated, with 27% in , and 2.4% in states with . 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 Koelemeij2012 (), we first ignore hyperfine structure and obtain transition dipole moments given by
(4) 
where denotes the radial wave function of nuclear motion obtained by numerical solution of the radial Schrödinger equation, stands for the internuclear separation, and denotes the permanent electric dipole moment function of the electronic ground state.
To take hyperfine structure into account we evaluate the dipole transition matrix element between two states and Brown2003 (); Koelemeij2012 (); Koelemeij2011 ()
(5)  
Equation (5) includes a transformation from the moleculefixed frame to the laboratoryfixed frame Brown2003 (), with denoting the polarization state of the electric field with respect to the laboratoryfixed frame ().
The linestrength is subsequently calculated by squaring and summing over the states Hilborn1982 ():
(6) 
with . The Einstein rate coefficients for spontaneous emission, absorption and stimulated emission are subsequently obtained as described in Tran2013 (), and used in Sec. 4.1 to calculate the expected (0,2)(8,3) spectrum.
3 Experiment
3.1 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 continuouswave (cw) 313 nm laser beam (see Koelemeij2012 (); Biesheuvel2013 (); Cozijn2013 () for details). The rf trap is placed inside an ultrahigh vacuum chamber with a background pressure of mbar. The rf trap operates at a frequency of 13.2 MHz, leading to Be radial trap frequencies of 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 electronimpact ionization as done by Blythe et al. Blythe2005 (), and monitored with a photomultiplier tube (PMT) and an electronmultiplied chargecoupleddevice (EMCCD) camera. EMCCD images show ellipsoidal mixedspecies 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 Koelemeij2012 ().
3.2 Spectroscopy of HD
The (0,2)(8,3) transition in HD is detected destructively through resonance enhanced multiphoton 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 copropagating 532 nm cw laser beam. Both lasers are directed along the trap axis and counterpropagate 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 =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 =2 are dissociated, leading to depletion of the =2 state (Fig. 2). During the remainder of the REMPD period, BBR repopulates the =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 frequencystabilized by locking its frequency to a nearby mode of a selfreferenced 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 shortterm stability, which is disciplined to the 1pps signal of a GPS receiver for longterm accuracy and traceability to the SI second with 210 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 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 socalled secular excitation Roth2006 (). 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 another 10s 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 transition are located in the range ( MHz, 140 MHz) around the hyperfineless frequency. As we expect Doppler broadening to 16 MHz, we divide this range into a set of 140 evenly spaced frequencies (spacing MHz) at which REMPD spectroscopy is performed. To convert possible timevarying 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 MHz, and using a few mW of laser power (corresponding to 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, , 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:
(7) 
Repeating this procedure for all 140 preselected 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 (with .
4 Results and Discussion
4.1 Spectral lineshape model
In order to determine the hyperfineless rovibratiotional (0,2)(8,3) transition frequency, , 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 , where is a fit parameter from which we later deduce the value of . Before fitting, the data are corrected for reactions with background gas (primariliy H). This procedure is described in Sec. 4.5.
We start with building a state vector which contains the population in all 62 hyperfine levels in the rotational states ranging from =0 to =5 with =0. This includes 97.6% of the total internal states of =0 given a BBR temperature of 300 K. We neglect the Zeeman splitting due to 0.19 mT Bfield 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 due to the Zeeman effect is considered separately in Sec. 4.6.2. Also the Stark effect, the electricquadrupole 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 =2 initial state interact with the 782 nm laser and BBR. We here make the simplifying assumption that any population in the 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 is obtained by solving the set of coupled rate equations:
(8) 
where and 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 offdiagonal elements describe the population a certain state gains from another. Since describes losses it only contains diagonal entries of the form
(9) 
where denotes the Einstein coefficient for absorption from a lower hyperfinerovibrational level to an upper level (which for the states and being restricted to those having and , respectively). The corresponding transition frequency is , is the intensity of the 782 nm laser, and represents a normalized response function averaged over the distribution of velocities of the HD ions. This involves an integration over all Doppler shifts observed by the HD ensemble, where is the wavevector of the 782 nm laser and the velocity in the direction. If the HD velocity distribution is thermal (Gaussian), depends only on the temperature in the direction, . However, the effects of micromotion and chemistry in the Coulomb crystal during REMPD lead to deviations from a thermal distribution. This implies that 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 are determined.
The matrix contains both diagonal and offdiagonal elements, taking into account the rate of exchange of population between all involved levels and mediated through all possible electricdipole transitions,
which correspond to spontaneous emission, stimulated emission and absorption by BBR, respectively. The relationship with Eq. (6) is given in Tran2013 (). denotes the BBR energy distribution function at frequency , which is given by:
(10) 
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 Koelemeij2007b ().
We use Mathematica to solve Eq. (8) in order to obtain . We subsequently find the relative HD loss, , by summing over the hyperfine state populations (62 hyperfine states in ) before (i.e. ) and after () REMPD and computing:
(11) 
Here and 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 , and . In what follows we furthermore assume that K, which is the average temperature in our experimental setup.
The question arises what the relation is between the signal defined in Eq. (7) and the fractional loss defined above. In previous work it was assumed that and are interchangeable Roth2006 (); Koelemeij2007a (); Schneider2010 (); Koelemeij2012 (). In Appendix B we study the dependence of the signal on the initial number of ions 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 by use of a slightly nonlinear function,
(12) 
where is defined as the effective Be temperature along the axis during the initial secular scan of the REMPD cycle (see Appendix B). This means we have to use a fivedimensional fit function
(13) 
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 on a sufficiently dense grid of values (encompassing 270 values between and 150 MHz), (20 values between 1 and 20 mK), , and (9 values between and Wm), which we interpolate (again using Mathematica) to find a fast, continuous and smooth approximation to the function , , which we insert into the function . The result is a fivedimensional (5D) fit function which is suitable for nonlinear least squares fitting. A 3D projection of assuming fixed values of , and is plotted in Fig. 4.
The reason for treating 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 40m. 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 as a free fit parameter which represents the average 782 nm laser intensity for all data points. Similarly, the variables and cannot be determined accurately a priori and are treated as free fit parameters as well.
4.2 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. Blythe2005 (), our mixedspecies ion ensembles contain not only Be and HD, but also BeH, BeD, HD, and HD. In this paragraph we focus on the latter two species which are created during the HD loading procedure through the exothermic reactions
(14) 
and
(15) 
After loading, the excess HD gas is pumped out of the vacuum chamber within a few seconds, after which the background vapor resumes its steadystate composition (predominantly H). 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 chargetomass ratios 1:4 and 1:5, corresponding to HD and HD, 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 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 (HD, HD). An average number of 43 HD ions is obtained together with 60 ions of heavier species. We cannot determine the relative abundance of HD and HD, but previous observations indicate a branching ratio between Eqs. (14) and (15) of 1:1 and, thus, equal abundances Oka1992 (); Pollard1991 (). 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 shottoshot 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 massselective secular excitation (Eq. (7)).
For the treatment of effects of chemistry on the lineshape in Sec. 4.4 below, we consider two scenarios:

Scenario a: ,

Scenario b: , ,
where scenario b reflects the one sigma upper variation (which well represents the width of the upper tail of the ion number distribution).
4.3 Effect of micromotion
The rf quadrupole field of the ion trap inevitably leads to micromotion of ions with nonzero displacement from the trap 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 Berkeland1998 (). 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 phasemodulation 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 . 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 , can be determined through fluorescence measurements of a trapped string of beryllium ions using a modified version of the photonrf field correlation technique Berkeland1998 (). 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, , is varied, the axial micromotion amplitude scales linearly with , while the radial one varies as . The latter behavior stems from the dependent confinement and the concomitant variation of the Be radial displacement with . Thus, measuring the micromotion amplitude for various values of allows separating the radial and axial contributions.
To determine the micromotion amplitude, we use a similar setup as described in Berkeland1998 (). Photons detected with the PMT are converted to electrical pulses and amplified by an amplifierdiscriminator, which generates a START pulse at time . Subsequently, a STOP pulse is generated at time at the first downward zero crossing of the rf signal. A timetoamplitude 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 STARTSTOP time delays, employing 1ns bins in a range 076 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 , where MB stands for MaxwellBoltzmann) is given by:
(16) 
where is the Be mass, the Be temperature, the 313 nm laser intensity, the saturation intensity of the 313 nm cooling transition in Be, MHz the detuning of the 313 nm laser light, the wavevector of the 313 nm laser, and and 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 term can be written as
(17) 
where is the amplitude of Be micromotion along the direction of the laser wavevector and is a time offset.
We extract by fitting Eq. (16) to the acquired fluorescence histogram. Repeating this procedure for various values of , a list of data points of the form 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 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 (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 in Eq. (9) with sidebands at frequencies with amplitudes . Here, denotes the wavevector of the 782 nm laser, and the are Bessel functions of the first kind, with an integer in the range .
4.4 Effects of chemistry in the Coulomb crystal
During REMPD, H molecules in the background gas can react with the ions in the Coulomb crystal. Such reactions can be divided into two classes: (i) elastic collisions, and (ii) inelastic collisions, during which a chemical reaction or charge exchange occurs, and chemical energy is converted into kinetic energy. In general, the kineticenergy transfer to the ion from elastic collisions with roomtemperature particles is much lower than the chemical energy released from inelastic collisions.
At close range , the electric field of the ion polarizes the neutral molecule which results in an attractive interaction potential . Here denotes the polarizability volume (in m)
If a neutral atom or molecule and an ion approach each other within a critical impact parameter , where and are the reduced mass and relative velocity of the pair, a socalled Langevin collision occurs, during which the collisional partners spiral towards each other and a chemical reaction can occur at very short range Hasted1964 (). 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 timeaveraged 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.
Number  Reaction  Energy of ionic product (eV)  Rate(s) 

1  HD + D + H  0.41  010 
2  HD + H HD + H  0.36  0.0042 
3  HD + H H + D  0.66  0.0014 
4  HD + H HD + H  0.016  0.0019 
5  HD + H D + H  0.017  0.0004 
6  HD + H HD + HD  0.022  0.0015 
7  Be(P) + H BeH + H  0.25  0.0019/0.005 
The rate (in s per molecule) of D production is dependent on REMPD time and frequency of the 782 nm laser. 
The rate is dependent on the fraction of time a Be ion spends in the excited P state, which is dependent on the 313 nm laser intensity and detuning ( MHz or MHz, respectively). 
The reaction numbered 1 corresponds to the REMPD process itself. From observations reported in Esry1999 () we infer that the ratio of HD D + H and HD H + D is approximately 1:1. The chargetomass 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 Pa in our apparatus Wineland1997 (). The different rate constants of reactions 2 and 3 illustrate the fact that HD can react with H in two ways: either the H 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 Oka1992 () and Pollard1991 () the probability of each scenario is approximately 50 %. In case the ion breaks apart, the probability that either a H or a D is donated to the H 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 in reaction 3 have the same chargetomass 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 Giles1992 (). 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 nonthermal. 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 velocity distribution shown in Fig. 9.
Details of the MD code are given in Appendix A. It turns out that the velocity distribution deviates from a Gaussian curve and is better described by a Gaussian Umarov2008 (), which is a Gaussian curve with higher wings parameterized by a parameter . For , the Gaussian is written as
(18) 
Here, and 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 , the function reduces to a regular Gaussian distribution. While we use Gaussians to describe the simulated data, the observation that 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 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 may shift by several hundreds of kHz by implementing a Gaussian with varying between 1.0 and 1.1. Note that another recent study of MD simulations independently confirmed the nonthermal character of velocity distributions of lasercooled ion crystals due to collisions with background gas molecules Rouse2015 ().
In the remainder of this section, we determine the values applicable to our REMPD measurement with the help of MD simulations. The value of is dependent on the number of trapped fast ions, . The more fast ions, the higher the value. We note that 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 Chen2013 (). 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 can be obtained by solving the rate equation:
(19) 
The term is the production rate from a number of particles, such as Be or HD. In Sec. 4.2 we determined =750 and or . is the rate at which fast ions are cooled and become embedded within the Coulomb crystal, while is the rate at which (fast) ions escape from the trap. The values of correspond to the reaction rates in Table 1. However, obtaining and 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. and 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 .

Scenario 2: A maximum number of fast ions is present in the trap, which is realized by setting and 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 .
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 HD and H occur less frequently in the trap, and in line with the extreme scenarios introduced above we assume zero HD and H for scenario 1, and 3 HD and 1 H 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 results.
For all possible combinations of fast ions present during REMPD (e.g. 1 D, 2 BeH and 1 HD or 3 D, 3 BeH and 3 HD ) an MD simulation is carried out. From these simulations, the HD velocity distribution is determined and a Gaussian is fitted which results in one 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 timedependent number of available HD ions in the target state. To take this properly into account we introduce a time and frequency dependent parameter as follows. For each of the four scenarios, the number of fast D is simulated on a grid of different REMPD durations, , and of different frequencies, of the 782 nm laser. On each point of this twodimensional grid, the number of fast D is combined with the number of other fast ions, and the corresponding value of is looked up in a library of values, obtained from many MD simulations performed with various combinations and abundances of fast ion species. Interpolation of the grid leads to a smooth continuous function , which is subsequently inserted into the lineshape function used in Eqs. (8) and (9). Figure 10 shows the 3D plots of for the different scenarios.
Besides the value, also the ion temperature 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 and share the same origin (namely collisions with fast ions), and in Fig. 11 we show the relation between and , obtained from fitting
Gaussians to simulated velocity distributions. We find a linear dependence of the form , with the slope 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 estimated temperature difference is 2.5(5) mK, where the uncertainty of 0.5 mK is treated as one standard deviation. The and dependent temperature is also included in .
A schematic overview of the various steps involved in the construction of the fit function , as well as the role of the four scenarios, is presented in Fig. 12.
4.5 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 molecules in the vacuum setup. In order to correct the spectroscopic signal for these socalled 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 which is detected as a small difference between the secular scan peak areas and . During each measurement session, a series of 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 is extracted from the relation
(20) 
with s and mapping onto (see Appendix B). The values of are inserted into a modified set of rate equations, which include the process of background loss reactions:
(21) 
Here, the vector equals of Eq. (8) extended by two additional rows which describe the occurence of ions in the form of HD or H. is a diagonal matrix containing which describes the HD losses. We take into account the fact that conversion of HD to H in reaction 3 in Table 1 is not detected by the method of detection through secular excitation because HD and H have the same masstocharge ratio. This also means that the measured values of 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 of 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 (see Sec. 4.1) is fitted to the corrected data points with free fit parameters , , and . To find an improved estimate of the background signal, the thus found values of , , 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 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 and . 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 chargetomass ratio 1:3 itself. For example, the conversion of HD to HD (reaction 2) and to H (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 and are at the level of a few percent, and are applied through Eqs. (11) and (20) (see also Fig. 12).
4.6 Spectrum, systematic effects and final result
The function (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 =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 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.
Scenario  

fit parameter  1a  1b  2a  2b 
(MHz)  0.46(33)  0.37(34)  0.03(32)  0.12(33) 
T(mK)  10.9(8)  11.0(8)  10.6(8)  10.35(80) 
I (10 Wm)  0.99(14)  1.04(15)  0.95(14)  0.98(15) 
(K)  2.79(3)  3.85(5)  2.82(5)  3.84(5) 
The correlation coefficients of the fit parameters are presented in Table 3 and the values for are graphically shown in Fig. 15.
The error bars represent the 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 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 by taking the mean of these two values, while interpreting the mean singlefit 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 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 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 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 , which is treated as a 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 frequencylocked 782 nm laser with the optical frequency comb is counted during REMPD. We use the beatnote 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 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 to each of the data sets gives 100 different values of . A histogram of the resulting distribution of values is shown in Fig. 16. From the histogram we find a mean offset 0.5 kHz from the frequency value 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.
T  I  
1  0.609  0.483  0.049  
T  0.609  1  0.690  0.038 
I  0.483  0.690  1  0.579 
0.049  0.038  0.579  1 
The 782 nm laser has a Gaussian lineshape with a width of 0.5 MHz. The convolution of this lineshape with the (Gaussian) Dopplerbroadened 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 with no significant effect on its value.
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 subcomponents of individual hyperfine lines. This is done by diagonalizing the sum of the hyperfine and Zeeman Hamiltonians (using the 0.19 mT field value used throughout the experiments), after which the eigenvectors and energy values are used to compute the stick spectrum of the (v,L):(0,2)(8,3) transition. This procedure is similar to that followed in Refs. Koelemeij2007a (); Bakalov2011 (), 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 , forming a new composite (and Zeemanshifted) lineshape function . 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 fit results based on versions with mT with a version with zero field. It is important to note that the linear polarization of the 782 nm laser is practically perpendicular to the 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 fit results, obtained with the above model for several / intensity ratios, and the fit result we found previously assuming zero magnetic field. For the 0.19 mT field used in our experiment, a small shift of 0.017(3) MHz is obtained, with the uncertainty due to the possible maximum polarization imbalance.
Origin  Shift  Uncertainty  

(MHz)  (MHz)  (ppb)  
Resolution (statistical fit error)  0  0.33  0.85 
Uncertainty value  0.25  0.23  0.61 
Uncertainty  0  0.080  0.21 
Ignoring populations =6  0  0.032  0.083 
variation in spectrum  0  0.028  0.072 
Doppler effect due to micromotion  0.055  0.020  0.052 
Frequency measurement  0  0.010  0.026 
BBR temperature  0  0.005  0.013 
Zeeman effect  0.0169  0.003  0.008 
Stark effect  0.0013  0.0001  0.0004 
Electricquadrupole shift  0  0.0001  0.0003 
2 order Doppler effect  0  0.000005  0.00001 
Uncertainty  0  0.000001  0.000003 
Total  0.0182  0.41  1.1 
This is a shift with respect to a scenario with a zero effect of or micromotion, and serves to illustrate the size of the effect. The shift itself, however, is absorbed in the value of . 
The value of these shifts is actually nonzero but negligibly small, and therefore ignored here. 
The ac Stark shifts due to the 782 nm, 532 nm and 313 nm lasers are 869 Hz, 452 Hz and 8 Hz respectively. These values represent the shift of the center of 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 Koelemeij2011 () 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 Bakalov2014 (). The secondorder 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 from scenario 1a based on the model containing rotational states 5 with an extended version containing 6 states we find a shift of of 28 kHz, which we treat as the uncertainty due to the neglect of population in . The rateequation model furthermore includes the BBR temperature, which we estimate to be 300 K with an uncertainty of about 5 K, caused by daytoday 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 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 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 (see Eq. (1)), which is estimated to be 50 kHz Bakalov2006 (). Comparing fits with values differing by 50 kHz we find that this has a negligibly small effect of 1 Hz on the result for .
Frequency of the (v,L):(0,2)(8,3) transition
We take the average of the fit values obtained from scenarios 1a and 2a, corrected by the systematic shifts described above, to find the value of the (v,L):(0,2)(8,3) transition frequency, MHz. 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, = 383,407,177.150(15) MHz Korobov2014b (). The two main uncertainty contributions are the statistical fit error of 0.33 MHz and the uncertainty in the factor scenario of 0.23 MHz.
5 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 subp.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 ( MHz) differs by only 0.6 p.p.b. from the theoretically predicted value ( 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 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 1p.p.b. level, and a test of molecular QED at the level of , 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 Salumbides2013 (); Salumbides2015 (); Biesheuvel2015 (). Third, the result presented here can be used to obtain a new value the protontoelectron mass ratio with a precision of 2.9 p.p.b. Biesheuvel2015 (), for the first time from molecular spectroscopy as proposed already four decades ago Wing1976 ().
Our analysis clearly demonstrates that the firstorder Doppler effect is responsible for the largest contribution to the uncertainty. In fact, removing the firstorder 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 twophoton spectroscopy. For example, in Refs. Tran2013 (); Karr2016 () an experiment was proposed in which the line in HD is addressed through a twophoton transition with nearly degenerate photons. Using counterpropagating laser beams with a narrow linewidth, the LambDicke regime may be reached in the present apparatus, such that firstorder Doppler broadening is entirely eliminated, while all other systematic effects could be controlled below the 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 protonelectron mass ratio, and even contribute to the determination of several other fundamental constants including the Rydberg constant, the deuteronelectron mass ratio, and the proton and deuteron charge radii Karr2016 ().
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 DutchFrench 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.
6 Appendices
Appendix A Molecular Dynamics simulations
MD simulations are implemented in Fortran code in order to realistically describe the dynamics of trapped and lasercooled ions in the presence of the timedependent 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, , the timedependent force from the trapping field field, , and an optional rf electric field, , which drives the secular motion:
(22) 
The radial and axial part of are given by
(23) 
and
(24) 
where is the secular angular frequency in the direction. The constants , and depend on the trap geometry, which are determined through a finiteelement analysis performed with the software package SIMION. The forces exerted on each ion are calculated, and trajectories are obtained using the Velocity Verlet method Verlet1967 (). Dopplercooling is included at the level of singlephoton 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 neutralion 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 ns. However, if for any of the trapped particles the condition is met, where is the distance between particle with velocity and the nearest particle , the time step is reduced by a factor of 10. Likewise, the step size is increased if the colliding particles separate again and .
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 highenergy ions. This allows for a larger integration time step (10 ns) and, thus, faster MD simulations.
Appendix B Derivation of the nonlinear fluorescence function
The relative HD loss during REMPD, , is related to the spectroscopic signal through the nonlinear function , which we derive here. The fluorescence yield during a secular scan depends on the Be temperature, , and is described by the scattering rate formula integrated over a MaxwellBolzmann velocity distribution. Neglecting micromotion effects, we take the scattering rate defined in Eq. (16). During a secular scan in the experiment we use the values = 2300 MHz and , and in what follows we drop these variables from the function argument of . While performing the secular scan, the temperature varies, which leads to a fluorescence peak as described by Eq. (16). The spectroscopic signal is the relative difference between the areas under the fluorescence peaks (see Eq. (7)). We may rewrite the area, , as
(25) 
where is a constant taking into account the collection and quantum efficiencies of the PMT or EMCCD imaging system, denotes the duration of the secular scan (10 s), and stands for the ‘effective’ value of during a secular scan. The effective temperature