Sample preparation
We first laser cool a hot strontium beam to 1 mK using a Zeeman slower and three-dimensional magneto-optical trap (MOT) on the 1S0→1P1 32 MHz transition at 461 nm. The atoms are further cooled to a few μK with a three-dimensional MOT operating on the 1S0→3P1 7.5 kHz intercombination transition at 689 nm. About 2,000 atoms are then loaded into a cavity-enhanced one-dimensional optical lattice at 813.4 nm. The cavity mirrors are placed outside the vacuum chamber and the lattice light, generated with an injection-locked Ti:Sapphire laser, is stabilized to the cavity using the Pound-Drever-Hall technique, using a double-passed acousto-optic modulation as a frequency actuator.
Stable laser
The 87Sr sample is probed on the 1S0→3P0 1 mHz clock transition with a 698-nm diode laser, which is stabilized to 26 mHz using a 40-cm Ultralow-Expansion glass (ULE) cavity14,21. The cavity enclosure features bipolar temperature control, a passive heat shield, a double-chambered vacuum, active vibration cancellation and acoustic shielding. The stabilized laser passes through an independent acousto-optic modulator (AOM) to steer the frequency of the clock laser light reaching the atoms.
Atomic servo
The offset of the clock laser frequency relative to the clock transition is determined with Rabi spectroscopy. In this work, measurements use Rabi pulse lengths from 160 ms to 4 s. The excited state population fraction after clock spectroscopy is measured by counting the number of 1S0 ground state atoms using 1S0→1P1 fluorescence, repumping the 3P0 excited state population to the ground state and again counting the number of ground state atoms. To lock the clock laser to the atoms, two excited state population measurements are performed on the clock transition (one on each side of the resonance centre). The difference between these measurements is used as an error signal, which is processed by a digital proportional-integral-derivative (PID) controller to steer the laser frequency onto the clock transition resonance.
Lock-in measurements with the atomic servo
Many systematic uncertainties are measured using a digital lock-in technique. In this scheme, an experimental parameter is set at one value, the clock transition is interrogated and the atomic servo computes a frequency correction22. The same procedure is then performed for a different value of the experimental parameter, using a second, independent atomic servo loop. As the experiment alternates between these two states, data are recorded and time stamped. Demodulation occurs in post processing. In all cases we seek the difference between the resonance centres measured by these control loops.
Density shift
The use of spin-polarized ultracold fermions suppresses s-wave interactions among our atoms; however, p-wave interactions that shift the clock transition can be significant at high precision. This density shift is proportional to the atomic density and insensitive to temperature (due to its p-wave nature and the one-dimensional lattice18). The density shift is greatly reduced compared with our previous generation Sr clock, due to the use of a cavity-enhanced optical lattice21. To measure this shift, we perform a lock-in measurement by modulating the atom number and looking for a frequency shift. Extrapolating this result to an operating atom number of 2,000 and trap depth of 71 Erec (where Erec is the lattice photon recoil energy), we reach a density shift of (−3.5±0.4) × 10−18 (Fig. 5).
The overlapping Allan deviation shows the density shift averaging down for 2,000 atoms and U0=71 Erec. The atom number was modulated between 2,400 and 12,000 atoms. The error bars represent the 1σ uncertainty in the overlapping Allan deviation estimator.
Lattice Stark shift
A lock-in measurement is performed for different lattice powers to study the intensity dependence of the lattice Stark shift. We determine this shift as a function of the optical trap depth at the location of the atoms, U0, which is proportional to the lattice intensity. The value of U0 is determined from the trap frequency along the lattice axis, which is measured using resolved sideband spectroscopy.
Changing U0 also modulates the trap volume, which creates a parasitic density shift that can mimic a lattice light shift. A Gaussian density profile predicts that the density shift scales similar to NU03/2. In our system, we experimentally verify this relation with negligible uncertainty. To cancel effects of the density shift on this measurement, we modulate the atom number according to the NU03/2 scaling such that there is common-mode density shift cancellation. To further ensure that the density shift is removed, in post processing we remove data with the largest atom number fluctuations until the average differential density shift is well below the final measurement precision.
U0 is stabilized with a laser intensity servo by monitoring the cavity lattice transmission. The lattice frequency is locked to a Yb fibre laser optical-frequency comb referenced to the NIST maser array. Varying the lattice intensity and frequency, we find the magic wavelength where the clock shift is not responsive to changes in U0.
Drifting background magnetic fields can cause the atom’s quantization axis to vary with respect to the clock laser polarization. This creates a drifting ac Stark shift. To solve this problem, we run a background magnetic field servo1 during the ac Stark shift measurement.
At some level, terms nonlinear in U0 (such as hyperpolarizability and M1-E2 shifts) will be required to precisely model the lattice intensity. To measure these small terms, ref. 43 relied on the ability to obtain lattices as deep as 103 Erec, to achieve a large lattice intensity modulation amplitude. However, the measurement could have been susceptible to technical issues such as a noisy tapered amplifier43 used to generate lattice light or parasitic density shift effects, which could be significant for such large changes in lattice trap25. To check whether our data supports terms nonlinear in U0 to model the lattice light shift, we use an F-test44 yielding F=0.17 for 22 degrees of freedom (corresponding to unbinned data). Therefore, within our measurement precision, our data only supports a linear model (Fig. 6). We note also that all our lattice Stark shift measurements are made near the clock operating condition, with each data point reaching the statistical uncertainty at the 1 × 10−17 level. Together, these points determine the Stark shift correction at the 1 × 10−18 level for the relevant condition of our clock.
For ΔU=U2–U1, the colour coding represents the values of U2 used for each point. Data that has the same value of ΔU are averaged to produce the points in Fig. 2a. The error bars represent the 1σ uncertainty in each point.
If we were to assume significant hyperpolarizability, we can use our data to infer a hyperpolarizability shift coefficient of (0.3±0.3) μHz Erec−2. This is consistent with the value reported in ref. 43. We could also use the hyperpolarizability coefficient of ref. 43 to correct our data, resulting in a minimal increase in our total uncertainty (from 2.1 × 10−18 to 2.4 × 10−18). However, as our statistical tests do not justify hyperpolarizability, only linear behaviour is assumed in our quoted ac Stark shift.
Temperature sensors
The in-vacuum temperature sensors, Heraeus thin-film PRTs, are mounted on the end of borosilicate glass tubes sealed to mini vacuum flanges. PRTs are a well-established technology for accurate thermometry and are ultrahigh-vacuum (UHV) compatible. The PRTs are pre-qualified by cycling their temperatures between an ice melting point (temperature stable to 1 mK) and 200 °C, and then choosing sensors that shifted <1 mK over four cycles. Four-wire phosphor-bronze connections to the sensors are soldered to electrical feedthroughs in the flanges. The sensor resistance is measured with a bridge circuit, comparing the PRTs with a 1-p.p.m. resistance standard. Resistance measurements are taken with forward and reversed excitation currents for data processing that removes thermocouple effects. Electrical error is quantified in Table 2 of the main text.
The mounting structures were installed in a test chamber and hand carried on a passenger flight to Gaithersburg, Maryland, for calibration at the NIST Sensor Technology Division. At NIST, the sensors were calibrated by comparing them with standard PRTs, traceable to NIST’s ITS-90 temperature scale and accurate to 1 mK, using a water comparison bath with 1 mK temperature stability45. The temperature uniformity in the isothermal region of the bath is within 1 mK. As thin-film PRT calibration shifts are quasi-random, mechanisms that could affect the calibrations would cause the two sensors to disagree. Agreement between the sensors throughout the shipping and installation process strongly suggests that no calibration shifts have occurred. Thin-film PRTs are generally robust against calibration shifts due to impacts.
We deal with immersion error by a two-stage process. First, the test chamber is filled with pure helium and the sensors are calibrated to the standard PRTs. Data are fit to the Callendar van Dusen equation, RHe=R0(1+AT+BT2), where R0, A and B are fit parameters. The helium acts as an exchange gas, enabling radial heat exchange along the glass stem and suppressing immersion error. Second, we measure the sensor resistance under vacuum, Rvacuum, as a function of Tflange–Tprimary. To quantify immersion error, we fit Rvacuum–RHe=C(Tflange–Tprimary)+Δ, where C and Δ are fit parameters. These two equations are used to obtain Tprimary as a function of Rvacuum and Tbase. Sensor self-heating is studied by varying the excitation current and extrapolating the results to zero current.
The sensors are installed in the clock chamber using a gas backflow. After installation, sensor baking at 150 °C means that 1.0 mK uncertainty, from thermal cycling, must be added to Table 2. One of the sensors can be translated inside the vacuum chamber with an edge-welded bellows. For clock operation, this sensor is positioned 2.5 cm from the atoms to prevent coating with strontium. The temperature difference between the atom location and 2.5 cm away is (1.45±0.03) mK, which is included in Table 2.
The sensor translation measurements and temperature measurements throughout the inside of the BBR shield confirm that temperature gradients are small, indicating a well-thermalized environment. Compared with previous efforts1, temperature gradients in the clock chamber are now smaller, because greater care was taken to minimize heat sources inside the BBR shield. To quantify the non-thermal heat shift, we model the geometry and emissivities of the vacuum chamber1. We find that our simulation is insensitive to changes in the emissivity values, and that the non-thermal heat correction is bounded below the 1 × 10−19 level for our level of temperature uniformity. The non-thermal correction has been included in the ‘Static BBR’ entry of Table 1 rather than listed in Table 2.
Decay measurement
After population is driven to the 3D1 state (Fig. 4a), 689 nm fluorescence from the 3D1→3P1→1S0 cascade is collected with a photomultiplier tube and then read out and time binned (using a 40-ns bin size) with an SR430 event counter. This photon counting setup provides 0.4 ns of timing uncertainty.
Our statistics have confirmed that the noise in this measurement is Poissonian. Simulating the measurement with the appropriate noise process shows that our fits should be given Poisson weighting to correctly obtain the fit uncertainty.
Other simulations show that the fit does not accrue an appreciable bias due to the specific pulse shape when we use pulses shorter than 300 ns or when we remove data when the pulse is on from the fit. To ensure that this fit bias is doubly suppressed, we take both approaches. We take 0.1 ns as a conservative bound on the remaining uncertainty.
We have calculated the correction due to the 3D1 hyperfine structure to be at the negligible 0.001% level. Therefore, we choose 0.1 ns as a comfortable upper bound on this effect.
We quantify systematic bias from stray distributed-feedback (DFB) laser light by switching off the AOM used to pulse this laser, while attempting to scan the 3P0→3D1 transition. We are able to observe this transition with stray light for exposure times of hundreds of milliseconds. By simulating the results of this scan, we can put a small 0.01 ns upper bound on stray laser light effects. We put the same bound on systematic bias from stray 2.6 μm radiation originating from the ambient heat in our lab.
We study the measured decay rate as a function of atom number to check for density dependence. We confirm that the decay rate is constant in density within our precision using an F-test, comparing a constant with a model linear in density. With a value of F=0.045 for the statistic (where there are 11 degrees of freedom), this test indicates no density dependence.
Dc Stark shift
A background dc electric field can arise from various sources, such as patch charges29 or electronics27. We have only measured a significant background dc Stark shift along one direction. This axis passes through the two largest viewports and the centre of the MOT coils.
To combat possible changes in the dc Stark shift, we actively suppress this shift with electrodes placed on the two large viewports. We measure ν+, the total dc Stark shift with the applied field in one direction, and ν-, the shift with the applied field flipped in direction. The background field is proportional to (ν+ −ν−), which is processed by a digital Proportional–Integrator servo. The servo applies a voltage to the electrodes to null the background field. The nonlinearity of the shift in electric field means that shift measurements average down rapidly when the background field is well cancelled. We measure a low 10−20 level shift with an uncertainty of (−0.1±1.1) × 10−19 in 20 min of averaging time.
Probe Stark shift
We perform this measurement by locking two independent atomic servos to 20 and 180 ms π-pulses. By keeping the pulse area, which is proportional to the square root of the probe intensity, fixed at a π-pulse, we can perform low-noise measurements of the probe Stark shift, which is linear in probe intensity. To resolve the shift well, we perform a large amplitude probe intensity modulation using a motor to move a neutral density filter in and out of the clock laser beam path. Control measurements confirm that this filter does not introduce systematic bias.
To prevent issues with many-body effects that might shift the clock transition frequency as a function of atom number, we study the probe Stark shift with a clock operation atom number of 2,000. Extrapolating this result to an operating clock pulse of 1 s, we observe a probe Stark shift of (−3.2±1.7) × 10−20.
First-order Zeeman shift
The first-order Zeeman shift is greatly suppressed by averaging locks to the two mF=±9/2 stretched states26. A residual first-order Zeeman shift could occur if there is appreciable magnetic field drift in between clock interrogations. We combat this by employing active background magnetic field cancellation1.
The difference between the mF=±9/2 stretched state frequency measurements is proportional to the background magnetic field. Drifts in this difference indicate a residual first-order Zeeman shift. Averaging down this difference, we measure a first-order Zeeman shift of (−1.6±2.0) × 10−19.
Second-order Zeeman shift
We measure the second-order Zeeman shift by monitoring the atomic frequency shift while modulating between high- and low-bias magnetic field values. We then extrapolate the observed frequency shift to operating conditions, using the fact that the shift is proportional to the bias field squared. The second-order Zeeman shift is measured as a function of the frequency difference between the mF=±9/2 stretched states, Δνstretch, which is proportional to the bias field magnitude. For clock operation, Δνstretch=300 Hz.
Background field drift can change the direction of the bias field, creating a time-varying lattice tensor ac Stark shift that would affect the measurement. To prevent this, we operate a background field cancellation servo. In addition, we reduce the sensitivity to drifts by aligning the field and the clock laser polarization. This is done by minimizing the amplitude of mF changing σ transitions. With this setup, we put a 10−20 level upper bound on systematic bias from field drift.
We measure the second-order Zeeman shift coefficient, the shift normalized by 
, to be (−5.82±0.07) × 10−16 kHz−2. This number is an atomic property and is independent of a particular measurement; thus, we average this result with four other determinations of this coefficient1,25,46,47. The final value for the shift at Δνstretch=300 Hz is (−51.7±0.3) × 10−18. We use a reduced 
inflated uncertainty to account for non-statistical variations between these data points.
Other shifts
Line pulling occurs when off-resonant spectroscopic features can slightly shift the clock transition frequency. This can be caused by imperfect spin polarization leaving population in mF states aside from ±9/2, clock laser ellipticity causing us to drive σ transitions, or clock transition sidebands that result from tunnelling between lattice sites. Calculations and data allow us to put a conservative upper bound on this effect at 1 × 10−19.
The first-order Doppler effect is not present in an optical lattice probed along the lattice axis, where the optical phase of the lattice and that of the clock probe lasers are referenced to a common mirror. A second-order Doppler shift is, in principle, present, but it is estimated to be at the 10−21 level. We put a comfortable 1 × 10−19 bound on this effect.
Collisions with the background gases in our UHV chamber can shift the clock transition frequency. At normal operating vacuum pressure, the background gas is largely hydrogen. We use the model of ref. 48 to put an upper bound on this effect of 6 × 10−19.
Steady-state error in the atomic servo could shift the measured clock transition frequency. We average lock data and find a servo offset of (−5±4) × 10−19.
Clock operation uses an AOM to scan the frequency and pulse the intensity of the clock laser. Phase transients occurring when this AOM pulses would appear as frequency shifts in clock measurements. We study the AOM phase transients by looking at the beat of the first AOM order with the 0th order on a digital phase detector. We also calibrated the phase transients of the detector itself. Drawing on the analysis of ref. 49, we infer an AOM phase chirp shift of (6±4) × 10−19.
Statistical methods
To calculate the shift of a given record, we perform a post processing demodulation of the data to extract a signal. The shift represents the mean of this signal. The statistical uncertainty is calculated from the s.e.m. If the reduced 
, the statistical uncertainty is inflated by 
. To remove the effects of residual laser drift, which is highly linear, from lock-in measurements, we use ‘three-point strings’. This analysis involves processing successive triplets of frequency measurements in linear combinations meant to cancel linear drift50.
