A terahertz vibrational molecular clock with systematic uncertainty at the $10^{-14}$ level

Neutral quantum absorbers in optical lattices have emerged as a leading platform for achieving clocks with exquisite spectroscopic resolution. However, the studies of these clocks and their systematic shifts have so far been limited to atoms. Here, we extend this architecture to an ensemble of diatomic molecules and experimentally realize an accurate lattice clock based on pure molecular vibration. We evaluate the leading systematics, including the characterization of nonlinear trap-induced light shifts, achieving a total systematic uncertainty of $4.6\times10^{-14}$. The absolute frequency of the vibrational splitting is measured to be 31 825 183 207 592.8(5.1) Hz, enabling the dissociation energy of our molecule to be determined with record accuracy. Our results represent an important milestone in molecular spectroscopy and THz-frequency standards, and may be generalized to other neutral molecular species with applications for fundamental physics, including tests of molecular quantum electrodynamics and the search for new interactions.

Here, we extend the lattice clock architecture to trapped neutral molecules and characterize the systematic frequency shifts in a Sr 2 vibrational lattice clock, achieving a total fractional systematic uncertainty of 4.6 × 10 −14 , comparable to the earliest realizations of optical atomic lattice clocks [58][59][60]. We measure the pure molecular vibrational frequency to 13 digits, establishing one of the most accurately known oscillator frequencies in the terahertz (THz) band thus far [38,61]. We leverage this to characterize the electronic ground potential of the strontium dimer-originating from the van der Waals bonding of two closed-shell atoms-by determining the dissociation energy of 88 Sr 2 with an accuracy surpassing the previous record for a diatomic molecule [54].
The results described here may be applied to a wide class of molecules (e.g., hydrogen isotopologues [62]), enabling the refinement of molecular quantum electrodynamics calculations, tests of fundamental laws, and potentially open new pathways for THz frequency metrology [63][64][65].

II. VIBRATIONAL CLOCK
The basic scheme of the molecular clock is as follows. We operate the clock on the pure vibrational transition between the weakest bound and most tightly bound irrotational states, (v = 62, J = 0) → (v = 0, J = 0), in the X 1 Σ + g ground potential of 88 Sr 2 . Here, v and J denote the vibrational and total angular momentum quantum numbers, respectively. This pair of clock states offers the largest possible pure vibrational splitting (or clock frequency, f clock ∼ 32 THz) in the ground state for our molecule, and can be used to obtain the molec-  The optical lattice (brown arrow) off-resonantly addresses an isolated rovibronic state in (1)1u to induce magic trapping conditions. (b) Experimental setup. The upleg master laser is stabilized to a reference cavity using the Pound-Drever-Hall (PDH) technique, and its phase coherence is transferred to the downleg laser via a frequency comb. The molecules are held in the 1D optical lattice. Co-propagating clock lasers are delivered to the molecules via an optical fiber with active fiber noise cancellation (FNC). The spectroscopic signal derives from absorption imaging of X(62, 0) photofragments at a slight grazing angle relative to the lattice. A Rb microwave standard acts as a flywheel oscillator, linking the molecular clock to GPS time for the absolute frequency measurement. Further information is given in the main text and Appendices A and B. (c) Two-photon Rabi oscillations between the clock states driven at the operational probe intensities (filled circles, experimental data averaged over 8 consecutive runs, error bars represent 1σ uncertainties; solid red line, analytical fit to an exponentially decaying sinusoid). We observe lines as narrow as 11(1) Hz (inset, green squares). For clock operation, we perform Rabi spectroscopy with a 30 ms π-pulse duration (indicated by the black arrow), resolving 30(2) Hz linewidths consistent with the expected Fourier limit (inset, black open circles). Each point in the inset is a single shot of the experiment, and solid lines are Lorentzian fits. ular dissociation energy (D 0 ) given the binding energy of X(62, 0). As a direct transition between J = 0 states is strictly forbidden, we drive the clock transition via a Raman process using two diode lasers detuned from the intermediate excited state (1)0 + u (v = 11, J = 1). This pathway through a deeply bound (1)0 + u state offers favorable Rabi frequencies, which facilitated stimulated Raman adiabatic passage (STIRAP) transfer between X(62, 0) to X(0, 0) in our preceding work [66]. By contrast, weakly bound (1)0 + u states near the intercombination line (which we utilized in Ref. [50]) are expected to have negligibly small transition strengths to X(0, 0) due to poor Franck-Condon overlap [67,68]. The relevant potentials are shown in Fig. 1(a).
The measurements take place in a retroreflected onedimensional (1D) optical lattice at ∼1005 nm. Trapped samples of ultracold molecules are created by photoassociating laser cooled strontium atoms at 2 µK to the (1)1 u (v = −1, J = 1) rovibronic state. This efficiently produces X(62, 0) ground state molecules thanks to the large transition strength [67]. Molecules formed in the undesired J = 2 excited rotational state are photodissociated, and the remaining atoms are wiped out of the trap with resonant 461 nm laser light. Our detection scheme relies on state-selective photodissociation of X(62, 0) followed by absorption imaging of the slow-moving atoms. As this destroys the molecular sample, the entire sequence is iterated to scan the clock transition. Appendix A contains an elaboration on the state preparation.
Raman clock spectroscopy is deeply in the Lamb-Dicke regime for co-propagating probes along the axial direction of the optical lattice (Lamb-Dicke parameter η LD < ∼ 0.02). The upleg (or pump) master clock laser at 378 THz (793 nm) is stabilized to a high finesse ultra-low expansion reference cavity with a measured drift rate of 30 mHz/s that is compensated using a linearly-ramped acousto-optic modulator. The phase coherence of the upleg is transferred to the teeth of an erbium-fiber-laserbased optical frequency comb by actuating on its repetition frequency [ Fig. 1(b)]. The downleg (or anti-Stokes) clock laser at 410 THz (731 nm) is phase locked to the comb, thereby inheriting the phase stability of the upleg. The carrier-envelope offset frequency of the comb is stabilized to a Rb standard that serves as the laboratory timebase. Since the Raman transition samples the correlated frequency difference of the clock lasers, spectral broadening due to the linewidth of the upleg master laser is greatly suppressed. The upleg is passed through an acousto-optic modulator (AOM1 in Fig. 1(b)), whose first order diffraction is used to iteratively step the difference frequency of the clock lasers across f clock . AOM1 controls the interrogation duration by pulsing the upleg, and we leave the downleg constantly irradiating (but blocked with a mechanical shutter during the state preparation process). Both clock lasers are delivered to the molecules via the same optical fiber, and active fiber noise cancellation [69,70] is implemented separately on each leg using the same phase reference surface; see Appendix B. Figure 1(c) shows two-photon Rabi oscillations driven by the clock lasers at the operational Rabi frequencies.
Using pulse durations of ∼100 ms, our apparatus is capable of producing clock lines with full width at halfmaximum as narrow as 11 (1) Hz, corresponding to a Qfactor of 2.9 × 10 12 (solid green squares in the inset to Fig. 1(c)). In the absence of other fields, the blackbody radiation (BBR) limited lifetimes of the clock states exceed 10 5 years in a room temperature environment [50], suggesting no fundamental limit for the Q-factor. Nevertheless, the main technical challenges in the current iteration of the molecular clock are two-body molecular losses close to the universal rate [66] and lattice light-induced one-body losses for X(0, 0) that scale quadratically with the trap depth [50,67]. At the operational density and lattice trap depth in this work, these losses quench the spectroscopic signal fast enough that the molecular densities vary significantly over pulse durations of > ∼ 60 ms. As  Fig. 1(c)). A typical spectrum consists of 15 experimental iterations (taking a total duration of ∼20 s), from which we determine the line center by fitting a Lorentzian function.

A. Systematic evaluation
Table I details the uncertainty budget of the molecular clock under the operational conditions of this work. Summing the uncertainties of all contributors in quadrature, we report a total systematic uncertainty of 4.6 × 10 −14 .
We leverage the short-term frequency stability of our reference cavity to average down the uncertainty of a given systematic. Most frequency corrections in Table I are determined by probing the clock transition in an interleaved fashion; i.e., we alternate an experimental parameter between two values and record the corresponding pair of line centers [71]. This is repeated to gather statistics, and the shift in the line center is computed as a weighted average. The clock shift, ∆f clock , is defined as the frequency shift relative to the unperturbed clock. The clock shifts are extrapolated to determine the frequency correction for the clock at the operational parameter value. We quote the revised (or scale-corrected) standard errors for all fit parameters obtained from weighted fitting or averaging; i.e., we scale up the statistical uncertainties of the fit parameters by the square root of the reduced chi-square statistic (χ 2 red ) if χ 2 red > 1, which indicates an overscattered dataset relative to the fitted model.

Lattice light shift
Magic-or state-insensitive-trapping conditions can be engineered for the vibrational clock states by offresonantly addressing X(0, 0) → (1)1 u (9, 1) with the lattice. In the present work, the protocol involves predominantly tuning the polarizability of X(0, 0) to match that of X(62, 0) at a magic detuning of 4.493(3) GHz [66]. Owing to a favorably large Franck-Condon factor, the transition strength between X(0, 0) → (1)1 u (9, 1) is the largest among all X(v, 0) → (1)1 u (v , 1) transitions [66,67], resulting in a magic detuning that is nearly 4× larger than that in our initial demonstration of the molecular magic wavelength technique [50], putting less stringent requirements on the frequency stability of the lattice laser.
Importantly, the neighboring (1)1 u (v , 1) rovibronic resonances are spaced at intervals of ∼2 THz, and may cause deleterious shifts due to lattice light impurity (e.g., amplified spontaneous emission (ASE) [72]). To mitigate this, the lattice light derives from a Ti:sapphire laser stabilized to the same optical frequency comb described in Section II. This also permits the lattice frequency, f latt = c/λ latt , to be determined with kHz-level accuracy. To further suppress ASE impurity at the magic detuning, the light is filtered through a linear cavity (finesse of 50, and free spectral range of 2.9 GHz) before delivery to the experiment by a single-mode polarization maintaining fiber. A stable, weak reflection from the vacuum window is used for lattice intensity stabilization during normal operation. The linear lattice polarization defines the quantization axis for the magnetically insensitive X 1 Σ + g states [73]. We investigate the effect of lattice light over a range of f latt . At each f latt we make interleaved measurements of the clock shifts, alternating the trap depth U 0 between a reference depth and four other depths spanning from 300 E r to 1100 E r , where E r ≡ h 2 /(2M λ 2 latt ) is the recoil energy and M is the molecular mass. The trap depths are determined from the axial trapping frequencies (Appendix C). Small corrections (< 0.3×10 −14 ×f clock ) were made to account for density shifts. As shown in Fig. 2(a), our measurements reveal nonlinear light shifts as a consequence of molecular hyperpolarizability. The higherorder transitions that account for this effect will be investigated in future work, but we hypothesize a connection with previously observed quadratic lattice scattering rates in a similar experiment [50].
In order to characterize the lattice light shifts, we adopt the thermal model described in Ref. [74] and write the clock shifts as where α * and β * are empirically obtained from parabolic fits to the measured differential shifts. These parameters are effective values dependent on the trapping conditions: α * is related to the differential electric-dipole (E1),  ) with respect to a reference trap depth (∼ 500 Er), and fit the data to parabolas (solid lines) with a global quadratic parameter, −β * . (b) Linear light shift coefficient, α * , versus lattice frequency (color code matches (a)), and the linear fit (black solid line). α * is predominantly due to the E1 differential polarizability and is nulled at fzero. By tuning α * , we can find conditions where the sensitivity of ∆f clock to fluctuations in U0 is minimal at our operational trap depth of 487(4) Er (dark green points). Error bars represent 1σ uncertainties. magnetic-dipole (M 1) and electric-quadrupole (E2) polarizabilities, while β * is related to the differential hyperpolarizability. Crucially, the polynomial form of Eq. (1) hinges on a linear scaling of the sample temperature with U 0 , which we verify to hold true for our molecules using Raman carrier thermometry (Appendix C). We do not expect non-polynomial terms [75] to be significant at the level of the current evaluation.
The fits give β * = −6.81(22) × 10 −5 Hz/E 2 r as a global parameter. Additionally, the results for α * versus f latt are shown in Fig. 2(b), and a linear fit yields a sensitivity slope ∂α * /∂f latt = −0.0796(16) Hz/(MHz E r ) as well Under these conditions, ∆f clock is first-order insensitive to changes in U 0 (dark green points in Fig. 2).

Probe light shift
Probe light shifts pose an inherent challenge for twophoton spectroscopy. This is even more so for scalar clock states (J = 0) that preclude the use of laser polarizationbased cancellation schemes [76]. Here, the clock shifts scale linearly as the probe intensities are low, and are related to the differential polarizability at the respective probe wavelength (λ p ), where α v is the E1 polarizability for the vibrational state v, I p is the probe laser intensity, and p ∈ {↑, ↓} specifies the laser: upleg (↑) or downleg (↓). Figure 3 shows that linear extrapolation of probe shifts suffices for a molecular clock at the few 10 −14 level. While tailored pulse sequences to alleviate probe light shifts have been proposed [77][78][79][80], for this evaluation we opted for a more straightforward strategy. We can minimize the total probe light shift by using socalled balanced intensity ratios satisfying the condition . At the same time, a large Raman detuning-relative to the intermediate (1)0 + u (11, 1) excited state-is preferred so that off-resonant scattering from the probes have a negligible effect on the accessible coherence times. Figure 3 demonstrates that such conditions exist in our clock for blue detunings where the baseline polarizability differences at the probe wavelengths have opposite signs, in agreement with our polarizability model (Appendix E). We operate at a Raman detuning of +14.973 GHz, much greater than the 5 MHz natural linewidth of the intermediate state [66]. We evaluate ∆f clock for each leg separately. Using a motorized neutral density filter, we switch between two intensity values for one leg while keeping that of the other leg at its operational value. The π-pulse durations are adjusted accordingly. Typical settings for the interleaved measurements are (P ↑,0 , 9P ↑,0 ), and (P ↓,0 , 3.5P ↓,0 ) where P p,0 = I p,0 (πw 2 p /2) are the operational powers measured with a calibrated power meter immediately before the vacuum window. These shifts are scaled by the measurement lever arms to obtain the clock corrections at the operational settings: −(∆f clock /∆P p )×P p,0 . We find the corrections to be −277.5(1.4) × 10 −14 for the upleg, and 309.0(1.7) × 10 −14 for the downleg. Drifts in ∆P p are at the sub-percent level over the ∼2000 s duration for each probe light shift evaluation, and the weighted averages of f clock typically have χ 2 red ∼ 1. Accurate knowledge of the beam waists w p is not necessary as they are robust during an evaluation, and they are common factors that drop out in calculations. Long-term drifts due to beam pointing instability may be monitored and countered by benchmarking the probe intensities using the molecules (e.g., using an Autler-Townes frequency splitting, an onresonance scattering rate, or the two-photon Rabi oscillation frequency), which we leave to future work.

Blackbody radiation shift
Homonuclear dimers are infrared inactive, conferring natural insensitivity to blackbody radiation (BBR). Using the formulas derived in Ref. [81], the frequency correction due to BBR is calculated to be −0.70 (14) Hz at an effective temperature of T eff = 303(5) K; see Appendix D for a description of the chamber thermometry. The uncertainty is dominated by ab intio calculations of the dc polarizabilities of the clock states (Appendix E). Comparison with experimentally measured ac polarizability ratios show agreement at the level of 10-20%, to be expected from typical accuracies of theoretical transition strengths. Therefore, we assign a conservative fractional uncertainty of 20% for the BBR shift.

Density shift
Due to their bosonic character, our 88 Sr 2 molecules are unprotected against s-wave collisions. The onedimensional lattice forms a series of microtraps, each with a trap volume proportional to (T /ω 2 ) 3/2 . Here, T is the temperature of the molecules, andω is the geometric mean of the angular trapping frequencies. We investigate density dependent shifts arising from dimer-dimer collisions by modulating the average number of molecules per lattice site (N mol/site ) at the beginning of the clock pulse. This is achieved by inserting a wait time immediately after photoassociation (PA) so that two-body collisions naturally reduce the molecule number [66,67]. Fluctuations in N mol/site are typically < 20%, and we assume equal occupancy across filled sites. Since both T andω 2 scale similarly with U 0 , and the lattice intensity is stabilized, N mol/site is a robust observable proportional to the molecular density.
Assuming linear density shifts, we scale our differential measurements to find ∆f clock at the normal operating value of N mol/site = 1. Figure 4(a) summarizes the measurements performed at various number differences (∆N mol/site ) suggesting a correction of −0.20(10) Hz, or −0.63(31) × 10 −14 in fractional units, due to collisional shifts. Control measurements using spectra taken under common experimental settings do not show evidence of spurious offsets in our data [ Fig. 4(b)]. It is instructive to compare the size of our density shift with similarly performing atomic clocks. From a trap calibration (Appendix C) we estimate that the shift coefficient has a magnitude of 2.9(1.5) × 10 −25 cm 3 after normalizing by the transition frequency. This is rather similar to the analogous optical atomic clock with bosonic 88 Sr (∼ 2 × 10 −25 cm 3 [82]), while being orders of magnitude smaller than in Cs (∼ 1 × 10 −21 cm 3 [83,84]) or Rb microwave clocks (∼ 5 × 10 −23 cm 3 [85]).

B. Absolute frequency evaluation
As illustrated in Fig. 1(b), we reference all RF frequency counters and direct digital frequency synthesizers (DDS) in the experiment to a free-running Rb microwave standard (our local timebase). Calibration of this Rb clock is accomplished by comparing its 1 pulseper-second (PPS) output with that of a dual-band global navigation satellite system (GNSS) receiver using a time interval counter (TIC); see Appendix G. The Rb clock serves as a flywheel oscillator to access Global Positioning System (GPS) time.
Each measurement trial of the absolute clock frequency is performed under operational conditions, where the molecular clock systematics are controlled at the level quoted in Table I. We repeatedly scan the clock transition to obtain a time series of the line centers, while simultaneously counting the repetition rate of the frequency comb. The probe light shifts were evaluated every trial to account for potential daily variations in probe laser beam pointing. The correction due to gravitational redshift is given in Appendix F. Figure 5 shows the results of the measurement campaign, consisting of 10 trials performed on separate days. A weighted average yields the absolute frequency of the 88 Sr 2 vibrational clock to be f clock = 31 825 183 207 592.8(5.1) Hz, with a fractional uncertainty of 1.6×10 −13 .

IV. CONCLUSION
Few frequency standards currently exist in the THz band [38,61]. Our molecular clock serves as a THz reference and can generate stable radiation at 9.4 µm via photomixing [86,87]. Alternatively, transitions in heteronuclear isotopologues could be driven directly with quantum cascade lasers [88,89]. To our knowledge, f clock represents one of the most accurately measured pure molecular vibrational frequencies to date. The fractional uncertainty is on par with that of the unidentified rovibrational interval in OsO 4 near the R(10) (00 0 1)-(10 0 0) emission line of the 12 C 16 O 2 laser. This absorption line in OsO 4 is a secondary representation of the SI second [61], and  Table I was compared directly against a primary cesium standard by stabilizing a CO 2 laser to the specific saturated absorption feature of OsO 4 in a high-finesse cavity [90,91]. We expect to reduce the uncertainty of our local timebase calibration to the same level as the molecular clock systematics (or better) by upgrading to a standard with intrinsically lower instability and utilizing two-way time transfer schemes.
Molecular spectroscopy is increasingly appreciated as a fertile ground in the search for new physics. The reported Hz-level molecular clock is a starting point for elucidating the bonding of the Sr 2 dimer across a large range of internuclear distances. The isotopologues of Sr 2 have different nucleon numbers, and comparison of their vibrational spectra may permit the investigation of hypothesized hadron-hadron interactions [32].
Gaining access to longer coherence times is a general strategy for improving the systematic uncertainty. It would enable the excitation of narrower molecular resonances, expediting the evaluation of a systematic shift. Operating at lower trap depths would considerably suppress the lattice light-induced one-body losses of the deeply bound vibrational state, X(0, 0). To this end, we plan to reorient the lattice in a future upgrade such that its tighter axial dimension is along gravity to permit the confinement of molecules in a shallower trap. Notably, atomic lattice clocks have entered the ∼ 10 E r regime [96], and adopting these recent techniques should further mitigate the lattice light-induced losses. This may be supplemented by deeper atomic cooling [97,98] prior to photoassociation. Moreover, given that the lattice light shift is the most significant systematic in this work, operation at shallower trap depths would directly improve the clock accuracy. Longer Rabi interrogation times imply smaller effective Rabi frequencies, thus the operational probe laser intensities can be reduced, which lowers the total probe light shift. Future work may circumvent collisional shifts and losses by preparing samples with single molecule occupancy in a three-dimensional lattice [99][100][101][102], or an optical tweezer array [103][104][105][106].
In summary, we have demonstrated a vibrational molecular clock with a total systematic uncertainty of 4.6 × 10 −14 , entering a new domain in high-resolution molecular spectroscopy. Our results are enabled by merging the key strengths of atomic clock techniques with molecular quantum science.

ACKNOWLEDGMENTS
We gratefully thank J. Sherman and T. H. Yoon for advice on the absolute frequency measurement and insightful discussions, and V. Lochab for early contributions to the vacuum chamber thermometry. We also thank the anonymous referees for their careful reading of the manuscript and useful suggestions. Our experiments start with a thermal beam of 88 Sr atoms decelerated by a Zeeman slower and laser cooled in a first-stage (blue) magneto-optical trap (MOT) using the 1 S 0 → 1 P 1 transition at 461 nm. We repump on the 3 P 2 → 3 S 1 and 3 P 0 → 3 S 1 transitions at 707 nm and 679 nm respectively. A second-stage (red) MOT on the narrower 1 S 0 → 3 P 1 intercombination at 689 nm further cools the atoms to a typical temperature of 2 µK. Throughout the cooling sequence (∼ 500 ms), a one-dimensional optical lattice at 1005 nm overlaps with the atom cloud, and atoms with kinetic energies lower than the trap depth are loaded into the lattice. The surrounding magnetic field is lowered to <0.6 G to prepare for molecule production and spectroscopy. The lattice is formed by retroreflecting the lattice laser beam and, in this study, oriented horizontally with respect to gravity due to practical limitations.
Spin statistics and molecular symmetry imply that only even values of total angular momentum J exist for ground state 88 Sr 2 molecules. While rotational factors favor the decay from 1 u (J = 1) to X(J = 0), a finite number of X(J = 2) molecules still form. Thus, to purify the gas, we photodissociate the X(62, 2) molecules 30 MHz above the 1 S 0 + 3 P 1 threshold, imparting more than sufficient kinetic energy to guarantee these photofragments leave the trap. We do this concurrently with photoassociation (PA) by adding a frequency sideband to the PA laser with an acousto-optic modulator. We use a PA pulse duration of ∼2 ms. The remaining atoms are blasted out of the trap with resonant 461 nm laser light. For the operational trap depth in this study, we prepare 6×10 3 molecules in the initial clock state X(62, 0), spread across approximately 520 lattice sites. To mitigate density shifts during clock operation, we further hold the molecules for a short duration (∼150 ms), leveraging on the natural two-body inelastic collisions to reduce lattice occupancy to 1 molecule per site, averaged over filled lattice sites. Finally, photodissociation of X(62, 0) near the 1 S 0 + 3 P 1 threshold produces slow-moving atoms that we absorption image to use as our spectroscopic signal [ Fig. 1(b)]. Mechanical shutters provide secondary shuttering of all laser beams (except the lattice) prior to entering the vacuum chamber, in addition to the fast primary shuttering performed using acousto-optic modulators. In this manner, we ensure that the molecules interact only with the Raman clock lasers and the lattice during clock interrogation. A complete account of our molecule production and detection methods can be found in Refs. [107,108].

Appendix B: Raman clock lasers
Our reference cavity is formed from two fused silica mirrors bonded to an ultra-low expansion glass spacer placed in a vacuum housing maintained at the measured zero-crossing temperature for the coefficient of thermal expansion. The cavity finesse is > 3×10 5 from ring-down measurements. The upleg diode laser serves as the master clock laser, and is stabilized to the cavity using the Pound-Drever Hall technique. We phase lock the repetition rate of an erbium-fiber-laser-based optical frequency comb directly to the upleg. Observations of the counted repetition rate against a Rb standard actively steered by a GPS disciplined oscillator for over a month prior to the campaign reveal a cavity drift rate of 30 mHz/s, which we compensate using an acousto-optic modulator in the optical path of the master laser to the cavity. The frequency synthesizer that performs this linear feedforward compensation updates every second. The residual linear drift of the master laser due to imperfect feedforward is approximately 3 mHz/s during the campaign, consistent with the observed drift of the molecular clock line centers over the same period after accounting for the comb teeth difference (3 mHz/s × [1 − (731 nm)/(793 nm)] ≈ 0.2 mHz/s). By phase locking the downleg diode laser directly to the comb, the comb acts as a transfer oscillator, and the phase fluctuations of the clock lasers are correlated. To suppress phase fluctuations due to the microwave synthesizers, the beats of the clock lasers with the comb are chosen to have the same sign and frequency; i.e., f b↑ = f b↓ in Fig. 1(b).
The clock lasers are injected into the same polarization-maintaining single-mode fiber and delivered to the adjacent optical table where the experiments take place. Since the wavelengths are sufficiently different that the laser beams may sample non-identical paths in a given refractive medium, active fiber noise cancellation (FNC) on each clock leg is implemented using independent phase actuators (acousto-optic modulators AOM2 and AOM3 in Fig. 1(b)). The voltage-controlled crystal oscillators (VCXO) provide the RF frequencies for AOM2 and AOM3, and the FNC beats are phase locked to the same RF reference derived from a direct digital frequency synthesizer (DDS). The phase reference surface at the experiment table is a single partially reflecting mirror, while the surfaces on the laser table are mounted on a common rigid pedestal post with the clock lasers approaching the surfaces in the same direction. To minimize the number of optical elements and unstabilized path lengths, the clock lasers interrogate the molecules from the opposite direction as the PA and photodissociation lasers. The total uncompensated path in air is approximately 50 cm. The polarizations of the probes are identical, linear, and parallel to the small applied magnetic field, but perpendicular to that of the lattice in this work. During the state preparation sequence described in Appendix A, both clock lasers are blocked by a mechanical shutter before the beams enter the vacuum chamber. The 1/e 2 beam waists of the upleg and downleg laser beams at the molecules are 89 (5) µm and 114 (20) µm, respectively.

Appendix C: Trap calibration and Raman carrier thermometry
The axial (f ax ) and radial (f rad ) trap frequencies for the molecules are measured with resolved-sideband spectroscopy at the operational magic lattice wavelength [ Figs. 6(a,b)]. To enhance the transition rates for the axial sidebands, we use counter-propagating probes to interrogate the naturally nearly-magic Raman transition between the adjacently bound vibrational states X(62, 0) → X(61, 0). To excite the radial sidebands, we use the Raman clock transition X(62, 0) → X(0, 0) but intentionally introduce a small relative misalignment in the probe beams. The trap depths are calculated using U 0 = M λ 2 f 2 ax /2. We note that, in general, U 0 /E r = M 2 λ 4 latt f 2 ax /h 2 , scaling as the square of the particle mass. Thus given the same trap frequencies and λ latt , U 0 /E r is numerically 4× larger for Sr 2 molecules than Sr atoms. Assuming a Boltzmann thermal distribution, the microtrap volumes are calculated as V = 2πk B T /(ω 2 M ) 3/2 , whereω ≡ 2π(f ax f 2 rad ) 1/3 and M is the molecular mass.
To determine T , the temperature of the molecular ensemble, we perform Raman carrier thermometry [66,109] with co-propagating probes scanning X(62, 0) → X(0, 0) at a non-magic wavelength (tuned >0.3 THz away from the operational magic wavelength to maximize the polarizability difference). As shown in Fig. 6(c), the dependence of T against U 0 is well described by a linear fit to the data.

Appendix D: Vacuum chamber thermometry
Short of finite-element modeling, we may make a basic estimate for the effective solid angle (Ω eff angle,i ) subtended by the i-th surface surrounding the molecular sample as Ω eff where η i is the surface emissivity and Ω angle,i is the geometric solid angle. The total effective solid angle is normalized to 4π. We use values for the emissivity of various materials from available literature [13,[110][111][112][113]. These are 0.91 for glass (fused silica), 0.54 for sapphire, and 0.08 for stainless steel. The sapphire window facing the Zeeman slower is heated to 430(10) K and subtends a geometric solid angle of 0.04 sr. Among the fused silica window viewports with direct line-of-sight to the molecules, there are 8 with a diameter of 33.78 mm, 4 with a diameter of 69.85 mm, and 6 with a diameter of 114.3 mm. The diameter of the spherical vacuum chamber is approximately 240(10) mm, and the surface area consisting of stainless steel is approximated as the spherical surface area subtracted by the total area encompassed by the viewports.
At the present level of precision, it is enough to estimate the temperature environment of the stainless steel vacuum chamber using four negative temperature coefficient (NTC) thermistors affixed to its exterior. The largest (smallest) sensor reading is T c,max (T c,min ). We model the temperature gradient as a rectangular distribution [114] and estimate the temperature of the vacuum chamber to be (T c,max + T c,min )/2 = 302 K, with an uncertainty of (T c,max −T c,min )/ √ 12 = 1 K. Conservatively, the fused silica windows are within ±2 K of the temperature of the stainless steel chamber. The line of sight from the molecules to the hot oven is blocked using an in-vacuum mechanical shutter during clock spectroscopy.
Following Ref. [115], the effective temperature (T eff ) that enters into the BBR shift calculation is such that where T i is the temperature of the i-th surface. We estimate an effective temperature of T eff = 303(5) K.

Appendix E: Clock state polarizabilities and BBR shifts
The electric dipole (E1) polarizabilities of vibrational states in the X 1 Σ + g potential are calculated using the sum-over-states approach [116]. For states with total angular momentum J = 0, the polarizability is a scalar quantity independent of the polarization of the electromagnetic field. Let |i represent the rovibronic wavefunction of X(v, J = 0) and |f represent the rovibronic wavefunction of a state that is E1-allowed from |i . The polarizability of a molecule in X(v, J = 0) is Here, d 0 is the component of the electric dipole moment operator parallel to the lab-frame quantization (Z) axis, ω f i is the angular frequency of the rovibronic transition, and ω is the angular frequency of the electromagnetic field. The dc (static) polarizability is recovered for ω = 0. The sum in Eq. (E1) is evaluated for bound-to-bound transitions to singlet ungerade excited potentials and bound-to-continuum transitions by discretizing the continuum. We include contributions from the 1 Σ + u potentials correlating to 1 S + 1 P and 1 S + 1 D, as well as the 1 Π u potentials correlating to 1 S + 1 P , 1 S + 1 D, 3 P + 3 P and 3 P + 3 D. The (1) 1 Σ + u potential is taken from the ab initio calculation in Ref. [117], while the doubly-excited (3) 1 Π u and (4) 1 Π u potentials were calculated using the multireference configuration interaction method (MRCI) with the MOLPRO package [118]. The remaining potentials (including X 1 Σ + g ) are empirical [95,119]. We omit spin-orbit and non-adiabatic couplings between the potentials. The convergence of our results is not changed by the inclusion of further high-lying potentials.
Laser wavelengths in the range 800-1200 nm can drive transitions from X 1 Σ + g to the short-range part of (1)1 u .
These singlet-triplet transitions are relatively weaker than singlet-singlet ones, but become important if the laser is tuned near a resonance; e.g., in the case of a magic wavelength optical trap. To properly account for these situations, we additionally include the Morse/Longrange potential of (1)1 u from Ref. [67] in the polarizability sum. Figure 7 shows the experimentally measured ac polarizability ratio α 0 /α 62 over a range of wavelengths, determined using a frequency-only method [66,67], along with the calculation using Eq. (E1) showing consistency within < 20%.
To calculate the BBR shift, we use the formulas derived in Ref. [81] and properly adapt them for molecular states. In Ref. [81], approximations were made to express the BBR shift in terms of the dc polarizability and a power series in k B T eff /(hω f i ). These approximations are also valid for our case since ω f i /(2πc) > 8000 cm −1 , corresponding to characteristic temperatures of >11500 K much greater than T eff . The so-called dynamic term contributes less than 0.5% to the total BBR shift. The correction to the vibrational clock frequency at T eff = 303(5) K is calculated to be −0.70 (14) Hz. Here, we quote a conservative fractional uncertainty of 20% based on the level of consistency between the measured and calculated ac polarizabilities [ Fig. 7].

Appendix F: Other clock systematics
Effects of magnetic field.-The use of singlet and irrotational X 1 Σ + g (J = 0) clock states confer a high degree of insensitivity to external magnetic fields. Hyperfine sublevels are absent in our dimer assembled from 88 Sr, which has a total nuclear spin of zero. While the excited molecular states of 88 Sr 2 near the intercombination line have been thoroughly studied and modeled in previous work [73,120], an equivalent quantum chemistry calculation of higher-order Zeeman shifts of X 1 Σ + g (J = 0) ground states is beyond the scope of this paper. Even without detailed theoretical modeling, we may experimentally investigate the extent to which our measurements are affected by magnetic field effects, including hypothetical Zeeman shifts of the clock states. We vary the applied magnetic field during clock interrogation with a lever arm of 3.2 G. The larger applied magnetic field slightly changes the photoassociation efficiency, which we partially compensate for by simultaneously varying the initial molecule number. Interleaved measurements obtain a differential shift of 0.05 (41) Hz. For hypothetical shifts that scale quadratically with magnetic field strength, the measurement suggests that these contribute < 5 × 10 −16 to the systematic uncertainty.
dc Stark shift.-The stainless-steel vacuum chamber is electrically grounded, and the molecules are held at a distance of approximately 120 mm from each of the fused silica viewports. We have operated the vacuum system for over a decade. Thus, we expect any stray charges on the viewports to have migrated and decayed to a negligible amount. Even if a hypothetical, improbably large voltage difference of 20 V were present between two opposite-facing viewports, using the dc polarizability difference computed with our model, the dc Stark shift multiplied by the number of such viewport pairs would amount to < 30 mHz, or < 10 −15 in fractional units.
Doppler shifts and phase chirps.-First-order Doppler shifts result from the relative motion of the lattice antinodes and the phases of the probe lasers. For example, this may be due to the mechanical motion of the lattice retroreflector, or phase chirps arising from the pulsing of an acousto-optic modulator (AOM) that diffracts a probe beam. Our upleg clock laser is pulsed by AOM1 before delivery to the molecules. A common solution in lattice clocks is to perform fiber noise cancellation of the probe(s) using the lattice retroreflector either as the phase reference surface or as a rigid support for a separate surface [121], which we will implement in future work. If uncompensated, AOM phase chirps can result in shifts as large as ∼100 mHz. We do not study phase chirps in this work. Consequently, we quote a conservative upper bound of 10 −14 for shifts originating from this effect. The second-order Doppler shift is < 10 −19 for the typical thermal speed of our molecule.
Lattice tunneling.-For the operational molecular trapping frequencies and temperature, we estimate that over 99% of molecules occupy motional quantum numbers n < 8, with 41% in the ground motional band (n = 0). A 1D lattice band structure calculation at the operational trap depth of U opt = 487 E r involving 1000 lattice sites (or equivalently, 1000 Fourier components) indicates that the bandwidth of n = 7 is 2 × 10 −5 E r , which translates to 0.02 Hz for our molecular mass (M ) and lattice wavelength (λ latt ). As a check, we verified that our calculation quantitatively reproduces the results of identical band structure calculations in available lattice clock literature [122,123]. We thus quote an upper bound on Doppler-like shifts due to the delocalization of the molecular wavefunction to be < 1 × 10 −15 .
Lattice light shift model.-We fit quadratic polynomials to the measured clock shifts for the lattice light shift evaluation to account for the observed hyperpolarizability light shifts. Given our trap frequencies, sample temperatures, and the linear scaling of sample temperature with trap depth, the polynomial fit is a reasonable approximation. The M 1-E2 shifts (that microscopically scale as √ U 0 [75]) are included in the α * effective parameter [11,13,74,124]. In future work, calculating the differential M 1 and E2 polarizabilities would help quantify the error associated with the thermal model.
To test if higher-order polynomial terms are statistically significant, we fit the lattice light shifts to a cubic polynomial, with the quadratic and cubic coefficients as global fit parameters. The data suggests a cubic coefficient of −1.2(9) × 10 −8 Hz/E 3 r . While the addition of a cubic term shifts the value of f zero , the estimated frequency correction, in this case, remains consistent with the applied correction in the main text within their un-certainties. As such, we limit our fitting to a quadratic polynomial for the present evaluation at the 10 −14 level.
Line pulling.-Under operational conditions, the radial trapping frequency is 311(2) Hz, which is 10× larger than the full width at half maximum of the clock resonances. The clock and the lattice laser beams are coaligned over several meters, and the radial sidebands are not visible during normal clock operation. We model the radial sidebands as two Lorentzian peaks centered at their expected detuning from the carrier, with amplitudes equal to the size of the typical shot-to-shot signal variation. To put an upper bound on the line pulling effect, we compare the difference in the carrier line center returned by fitting a typical clock spectrum with the sum of three Lorentzians (i.e., two radial sidebands and one carrier), versus the case using just a single Lorentzian (as in Fig. 1(c)). We estimate the line pulling error to be < 1 × 10 −15 .
Scan-and-fit error.-To estimate the effect of shortterm cavity flicker noise on our peak fitting, we fit a linear function to a typical time series of molecular clock lines totaling ∼ 3000 s (the typical duration for a single evaluation of a given systematic under interleaved clock operation). The magnitude of the linear coefficient is < 10 mHz/s. For the present experiment, it takes ∼20 s to scan out all 15 points that make up the clock spectrum. Therefore, we estimate an upper bound for the scan-and-fit errors to be < 6 × 10 −15 .
As mentioned in Appendix B, the months-long average linear drift of the molecular clock line due to imperfect feedforward compensation is approximately 0.2 mHz/s. The feedforward parameters were set beforehand and unchanged during the campaign. This longterm drift would contribute a systematic offset of magnitude 1 × 10 −16 , which is negligible for the current evaluation.
Gravitational redshift.-We determine the elevation of our apparatus to be 51(5) m above mean sea level, which corresponds to a redshift correction of −0.18 (2) Hz. This correction has been added to give the reported absolute frequency of the clock.  Table I, the Rb clock calibration, and the gravitational redshift (Appendix F). Error bars are not shown for visual clarity. Trial numbers on the top right corner of each plot correspond to those in Fig. 5(a). The histogram of all measurements is shown in Fig. 5(b). The axes labels are displayed in the plot for Trial 1 (top left), and are the same for the remaining plots.
comb, and the molecular clock systematic corrections).
f Rb /f GPS is the frequency of the Rb clock relative to GPS time. GPS time is closely steered toward UTC(USNO), and f GPS /f UTC(USNO) is the frequency of GPS time relative to UTC(USNO). Finally, f UTC(USNO) /f TAI is the frequency of UTC(USNO) relative to TAI, and f TAI /f SI is the frequency of TAI relative to the SI second.
The molecular clock is operated intermittently due to its complexity and the practical challenges of our present experimental apparatus. As such, the clock was not continuously phase-linked with the SI second. To address this, we expand Eq. (G1) as Here, T 1 is the typical up time of the molecular clock corresponding to one measurement trial [ Fig. 8], T 2 = quency is averaged. We have assumed that the SI second and molecular clock frequency are unchanging in time.
Over the length of the campaign, the scale intervals of GPS time, UTC(USNO), and International Atomic Time (TAI) differed by < ∼ 10 −15 , and the daily fractional changes in the frequency of GPS time relative to TAI are < ∼ 10 −14 [126]. Therefore, for the present study, except for f Rb (T 1 )/f GPS (T 1 ) and f clock /f Rb (T 1 ), all other ratios contribute a negligible uncertainty and may be assumed to be unity. This includes the extrapolation ratio, f GPS (T 1 )/f GPS (T 2 ), which is the frequency of GPS time during the molecular clock up times versus the frequency of GPS time broadcasted by the constellation over a month.
Each TIC measurement is started by the rising edge of the 1 PPS from the Rb clock and stopped by the rising edge of the 1 PPS from the GNSS receiver. Thus, the instantaneous fractional frequency offset of the Rb clock relative to the frequency of GPS time, r = [f Rb /f GPS ]−1, is quantified by the instantaneous slope of the logged TIC measurements as a function of elapsed time. This measurement is susceptible to noise in the satellite link, as well as the instabilities of GPS time and the Rb clock. Comparisons with an identical, independent free-running Rb clock indicate that the Rb clock reaches an instability flicker floor of approximately 3×10 −13 ≡ σ Rb after ∼ 5× 10 3 s of averaging time (comparable to typical durations of T 1 ), but worsens to ∼ 10 −12 for time periods over 24 hours. This poses a conundrum, because at least 24 hours of continuous averaging is typically required to achieve an inaccuracy and instability of < 10 −13 using one-way GPS time transfer [127], but the Rb clock is not a good flywheel on these time scales.
Rubidium microwave standards are more susceptible to unpredictable environmental perturbations than, for instance, hydrogen masers [128], making it challenging to construct a reliable noise model. Therefore, for every trial, we operationally extract r = [f Rb (T 1 )/f GPS (T 1 )]−1 through linear fitting of the TIC measurements as a function of elapsed time, restricting the fits to the durations coinciding with the up time segments of the molecular clock.
We judged a detailed characterization of the satellite link to be beyond the scope of this work. Geometric multipath effects and the diurnal variation in the ionosphere may introduce a systematic offset (≈ 2 × 10 −13 ≡ σ GPS,sys ), since a majority of the molecular clock up times were in the afternoon.
We estimate the fractional uncertainty of the extracted values of r to be σ 2 Rb + σ 2 GPS,tot , where σ 2 GPS,tot = σ 2 GPS,stat + σ 2 GPS,sys and σ GPS,stat ≈ 10 −13 × 86400/T 1 [s]. The uncertainty from linear fitting is an order of magnitude smaller. Occasionally, we manually realigned the Rb clock frequency relative to GPS if it exceeded a fractional offset of 1 × 10 −11 . This is not done during the molecular clock up times, nor within 24 hours of those segments to let the Rb clock settle. For every trial, we add a unique frequency correction r × [f GPS (T 1 )/f GPS (T 2 )] × [f clock /f Rb (T 1 )] to f clock /f Rb (T 1 ), obtaining the absolute frequency of the molecular clock.