# A lab-based test of the gravitational redshift with a miniature clock
network

Xin Zheng1,† 1University of Wisconsin-Madison, WI, USA Jonathan Dolde1,†
1University of Wisconsin-Madison, WI, USA Hong Ming Lim1 1University of
Wisconsin-Madison, WI, USA Shimon Kolkowitz1,∗ 1University of Wisconsin-
Madison, WI, USA

Einstein’s theory of general relativity predicts that when two clocks are
compared using light, a clock at a higher gravitational potential will tick
faster than an otherwise identical clock at a lower potential, an effect known
as the gravitational redshift [1]. This prediction has been tested by
comparing separated atomic clocks over a wide range of height differences [2,
3, 4, 5, 6, 7]. Here we perform a laboratory-based, blinded test of the
gravitational redshift using differential clock comparisons within an evenly
spaced array of 5 atomic ensembles spanning a height difference of 1 cm. We
measure a fractional frequency gradient of $[-12.4\pm 0.7_{\rm{(stat)}}\pm
2.4_{\rm{(sys)}}]\times 10^{-19}/$cm, consistent with the expected
gravitational redshift gradient of $-10.9\times 10^{-19}/$cm. Our measurements
can also be viewed as a demonstration of relativistic geodesy with millimeter
scale resolution. These results illustrate how local-oscillator-independent
differential clock comparisons can be harnessed for geodesy [9, 11, 10, 12,
8], and highlight their potential for emerging applications of optical atomic
clocks such as searches for new physics [13, 14, 15, 16, 17, 18],
gravitational wave detection [19], the realization of entanglement enhanced
clocks [20, 21, 22], and in explorations of the interplay between quantum
mechanics and gravity [23, 24].

Einstein’s theory of general relativity has thus far been consistent with
every experimental test performed [25], including classical [26, 27], modern
[28, 29] and strong field cosmological tests [30, 31, 32]. However, despite
the successful integration of special relativity and quantum mechanics as
quantum field theory, there is currently no theory that successfully unifies
general relativity with quantum mechanics. This motivates continued
experimental tests at new length scales, and suggests that performing
precision tests of general relativity with quantum systems may offer a way to
explore the interplay between general relativity and quantum mechanics [23,
24].

The gravitational redshift is a central prediction of general relativity.
Thanks to rapid advancements in their stability and accuracy [33, 34, 8, 35,
36, 37, 38], atomic clocks have now enabled tests of the gravitational
redshift over a wide range of length scales. Recent tests of the gravitational
redshift include a frequency comparison between two single-ion clocks with one
of the clocks elevated by 30 cm [3], comparisons between terrestrial clocks
and microwave atomic clocks in eccentric orbits which have produced the
strongest limits on deviations from the expected redshift [4, 5], a frequency
comparison between two synchronously linked portable 87Sr optical lattice
clocks with a 450 m height difference at the Tokyo Skytree tower [6],
resulting in the most precise terrestrial constraint on deviations from the
gravitational redshift at the $10^{-5}$ level, and a recent in-situ
synchronous frequency gradient measurement across a millimeter-scale atomic
ensemble with an unprecedented differential precision of $7.6\times 10^{-21}$
[7].

For two otherwise identical clocks experiencing the same gravitational field
with a height difference $\Delta h$, their frequency difference $\delta f$ due
to the gravitational redshift is given by

$\frac{\delta f}{f}\approx\frac{g\Delta h}{c^{2}}$ (1)

where $f$ is the clock frequency, $c$ is the speed of light, and $g$ is the
gravitational acceleration. Near the surface of Earth, this amounts to a
fractional frequency shift of $1.1\times 10^{-18}$ per centimeter of vertical
displacement. With optical clocks now reaching instabilities and inaccuracies
at the level of $10^{-18}$ and below [8, 35, 36], they are becoming a
sensitive probe of the point to point geopotential at the sub-centimeter
scale, where they are expected to outperform other methods of geodesy that do
not make use of spatial averaging [9, 11, 10, 12, 8]. For example, a blinded
comparison between two independent Yb optical lattice clocks was recently
performed with accuracy, instability and reproducibility all at the level
required to resolve sub-cm height differences [8]. In addition, a measurement
of the frequency gradient due to the gravitational redshift across a single
millimeter-scale atom ensemble was recently performed [7]. However, to our
knowledge, relativistic measurements of sub-centimeter scale height
differences between separate atom ensembles have not yet been demonstrated.

Emerging clock applications such as relativistic geodesy require transportable
optical clocks, which currently have poorer stabilities and accuracies than
those of laboratory-based clocks [39, 12, 11]. State-of-the-art laboratory-
based optical clocks often make use of bulky and immobile reference cavities
with second-scale coherence times [40, 41, 42] in order to achieve lower
levels of clock instability and to aid in rapid systematic evaluation,
limiting deployment in the field. It has recently been demonstrated that
differential measurements between single ions [43] and neutral atom ensembles
[44, 45], as well as differential spectroscopy by phase-coherently linking a
zero-dead-time optical lattice clock and a single ion clock [46], allow
interrogation times beyond the limit set by the local oscillators, opening up
the prospects for future applications with transportable or space-based clocks
[13, 19, 14, 15, 16].

Testing the gravitational redshift with a miniature clock network

In this work, we perform a local-oscillator independent, blinded test of the
gravitational redshift at the sub-centimeter scale using a spatially
multiplexed optical lattice clock [45] consisting of an array of 87Sr atom
ensembles trapped in a vertical, one-dimensional (1D) optical lattice (Fig.
1a). We prepare 5 atomic ensembles equally spaced by 0.25 cm, spanning a total
height difference of 1 cm. Synchronous differential comparisons are performed
between the 5 ensembles, resulting in 10 unique pairwise clock comparisons
recorded simultaneously, including 4 pairs at 0.25 cm, 3 pairs at 0.50 cm, 2
pairs at 0.75 cm, and 1 pair at 1.00 cm height difference. The gravitational
redshift is tested by mapping out the frequency differences between each
ensemble pair as a function of height difference.

The principles and basic operation of the multiplexed optical lattice clock
were described and experimentally demonstrated in Ref. [45]. However, several
improvements to our apparatus and experimental procedure have been made since
this prior work. A deeper initial optical lattice ($130\leavevmode\nobreak\
E_{\rm rec}$ trap depth, where $E_{\rm rec}$ is the recoil energy of a lattice
photon) during atomic ensemble loading and in-lattice cooling provide larger
atom numbers ($>2000$ atoms per ensemble) with reduced atomic temperatures
both axially ($>99\%$ occupancy in the lattice ground band) and radially
($<200$ nK). Synchronized Ramsey measurements are subsequently performed at a
shallower operational lattice depth ($u_{op}$) of $15\leavevmode\nobreak\
E_{\rm rec}$. Combined with common-mode rejection of clock laser noise, we
achieve 32 s atom-atom coherence times, more than 300 times longer than our
measured atom-laser coherence time of roughly 100 ms (Methods). This enables
differential instabilities below $1\times 10^{-17}/\sqrt{\tau}$ for all the
ensemble pairs, a factor of 3 reduction in instability over our previous work
($3\times 10^{-17}/\sqrt{\tau}$ with 6 ensembles). In addition, we now
suppress the residual magnetic field gradient along the lattice axis by a
factor of more than 10, reducing systematic uncertainties arising from Zeeman
shifts. These advances allow us to rapidly evaluate the differential
systematic shifts due to environmental perturbations and thus perform a
precision test of the gravitational redshift.

For our measurements, we utilize synchronous Ramsey spectroscopy in
conjunction with spatially resolved fluorescence imaging to probe the clock
transition along the ensemble array. The optical lattice is operated near the
magic wavelength where the differential polarizability between the ground
(${}^{1}S_{0},g$) and clock (${}^{3}P_{0},e$) states is zero [47, 48]. We
probe with an interleaved sequence using the magnetically insensitive
$\ket{g,m_{F}=\pm 5/2}\leftrightarrow\ket{e,m_{F}=\pm 3/2}$ clock transitions
[37], where $m_{F}$ is the projection of total angular momentum along the
quantization axis defined by the bias magnetic field. Taking the average of
the clock transitions with opposite sign nuclear spin states rejects first-
order Zeeman shifts and vector AC Stark shifts. The differential phase
$(\phi_{d})$ for each ensemble pair is extracted through least squares ellipse
fitting, and is related to the differential frequency $(\delta f)$ for the
pair through $\phi_{d}=2\pi\delta fT_{R}$, where $T_{R}$ is the Ramsey free
evolution time. A representative outcome from clock interrogation on the
$\ket{g,m_{F}=+5/2}\leftrightarrow\ket{e,m_{F}=+3/2}$ transition in shown in
Fig. 1b, where $(\delta f)$ is dominated by the first-order Zeeman shift
(Methods).

The differential frequency between each atomic ensemble pair includes a
contribution from the gravitational redshift as well as other frequency shifts
arising from differences between the two ensembles and their environments,
necessitating an evaluation of potential sources of systematic effects. To
avoid possible bias towards the expected result, we adopt a blinded
measurement protocol. A blinded offset gradient was randomly drawn from a
range of $\pm 5\times 10^{-18}/$cm, roughly 10 times the size of the expected
redshift gradient, and is automatically added to the extracted differential
phase by our data analysis code following the ellipse fitting step. This
blinded offset was not known to the authors during systematic evaluation and
data taking. The blinded value of the measured frequency gradient across the
array, the corrections for systematic shifts, and the systematic and
statistical uncertainties were all finalized prior to the removal of the
blind.

Systematic effects and error budget

The results of the full systematic evaluation are listed in Table 1. The
procedure for measuring the contribution of each potential systematic is
discussed in the Methods. Several effects dominate the extracted differential
frequency and the corresponding systematic uncertainty, and we highlight them
here. First, atomic interactions due to $p$-wave collisions between on-site
atoms lead to a frequency shift that scales linearly with atomic density [49,
50, 51]. In our differential measurement, the density shift is suppressed by a
factor of roughly 10 compared to the absolute shift by balancing the number of
atoms loaded into each ensemble. By intentionally varying the atom number
differences, the differential density shift at $u_{op}$ is evaluated to be
$-0.7(1)\times 10^{-19}$ per 100 atom number difference (Fig. 2a). In each
measurement run, the corresponding differential density shifts are corrected
for each ensemble pair individually.

The largest systematic shift in our system is the second-order Zeeman shift
due to the background magnetic field gradient ($\sim 1.5$ mG$/$cm). The
splitting of the transitions with opposite sign nuclear spin states provides a
measurement of the magnetic field difference between each ensemble pair, while
the overall magnitude of the applied bias magnetic field ($\sim 5.5$ G) is
measured independently using a more magnetically sensitive transition. This
bounds the uncertainty from the second-order Zeeman gradient to below $9\times
10^{-20}/$cm, limited by uncertainty in the second order Zeeman coefficients
for the clock transition (Fig. 2b and Methods).

The frequency shift due to black body radiation (BBR) is the dominant source
of systematic uncertainty for many room-temperature optical clocks[34, 8, 35,
36]. In our system, uncertainty arises due to temperature gradients in the
surrounding environment. Because the ensembles are arrayed vertically, the
primary contribution comes from differences in temperature between the top and
bottom recessed viewports of the science chamber, which are the closest
surfaces to the atoms. To evaluate the effect of this BBR gradient, we
intentionally introduce a thermal gradient by varying the temperature
difference between the two viewports by up to $\pm 1$ K. Mapping out the
resulting differential frequency shifts of the 10 ensemble pairs, we measure a
BBR gradient sensitivity of $-4.2(1)\times 10^{-18}/$cm per 1 K difference
(Fig. 2c and Methods). This results in a weighted average BBR gradient of
$-15.7(1.5)\times 10^{-19}/$cm under normal operating conditions, which is
determined from the temperature differences recorded with calibrated
temperature sensors that were continuously monitored during data taking. This
measurement also highlights the application of differential comparisons
between spatially multiplexed ensembles for BBR gradient evaluation at the
$10^{-19}/$cm level, which will be important for pushing the accuracy of room
temperature optical lattice clocks below the $10^{-18}$ level.

AC Stark shifts from lattice light contribute to the fractional frequency
uncertainty of state-of-the-art optical lattice clocks at the low $10^{-18}$
level [47, 48]. For differential comparisons between ensembles within a single
shared optical lattice, this is significantly suppressed. A differential
lattice AC Stark shift between two ensembles is caused by the relative trap
depth difference $\delta u$ arising from the lattice beam profile, and scales
linearly with $\delta u$ to first-order (Methods). In conjunction with the
multiplexed ensemble technique, we can rapidly map out both $\delta u$ and the
differential light shifts. In doing so, we also observe a residual spatial
light shift gradient of $-8(1)\times 10^{-20}/E_{\rm rec}/$cm, which depends
linearly on the lattice depth and the spatial separation between ensembles,
and does not depend on $\delta u$. We believe this gradient is likely due to
the differential tensor Stark shift arising from slight variations in the
orientation of the magnetic field vector across the ensemble array. This is
supported by the observation of a differential vector Stark gradient of
$-2.5(2)\times 10^{-18}/E_{\rm rec}/$cm in our system (Methods). Regardless of
the exact origin of the spatial light shift gradient, we are able to measure
and account for it by varying the lattice depth. Subtracting off this residual
spatial gradient, we observe correlations between $\delta u$ and the remaining
shifts (Fig. 2d), as expected. This allows us to extract the operational
lattice detuning from the effective magic wavelength [48, 7], and to
independently evaluate the lattice light shift gradient upon removal of the
$\delta u$ shifts (Fig. 2e).

Results and discussion

We performed 14 blinded measurements of the gravitational redshift under
normal operating conditions over a 3-week data taking campaign. In each
measurement run (ranging in duration from 1 to 4 hours), the frequency
gradient is determined by fitting a linear slope to the 10 measured
differential frequencies as a function of the pairwise height differences,
after taking into account systematic corrections such as ellipse fitting bias
and density shift. Corrections for spatially varying systematics such as the
lattice light shift, BBR shift, and second-order Zeeman shift are applied to
the measured gradient. Upon removal of the blinded offset gradient, we find a
weighted mean fractional frequency gradient of $[-12.4\pm 0.7_{(\rm stat)}\pm
2.4_{(\rm sys)}]\times 10^{-19}/$cm (where ‘stat’ and ‘sys’ indicate
statistical and systematic uncertainties, respectively), consistent with the
expected redshift gradient of $-10.9\times 10^{-19}/$cm within $1\sigma$ total
uncertainty (Fig. 3a). We note that the statistical uncertainty is scaled up
by the square root of the reduced $\chi^{2}$ statistic, $\chi^{2}_{\rm
red}=1.62$. Our measurement is consistent with a nonzero gravitational
redshift on the surface of the Earth for mm-cm scale height differences at a
confidence level of $4.4\sigma$.

Deviations from the gravitational redshift predicted by general relativity can
be parameterized by defining a modification parameter $\alpha$

$\frac{\delta f}{f}=(1+\alpha)\frac{g\Delta h}{c^{2}}$ (2)

to first-order of the gravitational potential difference [6]. We constrain
deviations from the predicted scaling by $\alpha=0.13\pm 0.23$ for millimeter
to centimeter scale height differences. We note that while the most stringent
constraints on $\alpha$ are at the $10^{-5}$ level [4, 5], those measurements
were performed at very different length scales, with height differences
roughly a factor of $10^{9}$ times larger than those used here. In addition,
we see a clear path to reducing the systematic uncertainty of our differential
comparisons by more than one order of magnitude (as was recently demonstrated
within a single ensemble [7]), and the techniques used here could be applied
in a larger apparatus or between ensembles in two separate vacuum chambers to
achieve spatial separations one the scale of 1-10 meters, opening the door to
laboratory-based tests of the redshift that are competitive with space-based
tests [4, 5] or tests with portable clocks [6].

As an alternative to the above analysis, we re-analyze the same experimental
data by taking the weighted average of the differential frequencies of each
ensemble pair from the 14 measurement runs, with systematic corrections
applied individually to each pair. The weighted mean differential frequencies
of ensemble pairs with the same height differences (0.25, 0.50, 0.75 and 1.00
cm) are shown in Fig. 3b. Through a final linear fit, we find a frequency
gradient of $(-11.9\pm 2.5)\times 10^{-19}/$cm, again fully consistent with
the expected redshift. Extrapolating the spatially varying systematic
uncertainties to 0 cm of separation, we find our gravitational redshift
resolution to be 1.3 mm, dominated by systematic uncertainty due to the
differential density shift, which could potentially be further reduced by
incorporating the recent technique for density shift cancellation [51] in
future work (Methods).

The measurement of gravitational potential differences between clocks with
sub-centimeter resolution is a major goal of relativistic geodesy [9, 11, 10,
12, 8]. In the preceding discussion and analysis we independently measured the
height differences between each ensemble pair and the local gravitational
acceleration $g$ and treated these values as known, but did not _a priori_
trust the gravitational redshift predicted by general relativity. From another
perspective, our measurements can be viewed as a proof-of-principle
demonstration of relativistic geodesy with millimeter scale resolution. Taking
as a given that the redshift is given by Eq. 1 and treating the ensemble array
as a network of spatially distributed clocks with unknown height differences,
we can extract the height ordering and relative height differences from the
measured gravitational redshifts and $g$, which was measured independently. We
find that we correctly assign the height order of the network, and that all of
the extracted height differences are within one standard deviation of the
known values (Fig. 4). This highlights the potential of synchronous
differential optical clock comparisons for geopotential measurements beyond
the geodetic limit [9, 11, 10, 12, 8], without placing stringent requirements
on the coherence time of the local oscillator.

Conclusion

We perform a blinded, precision test of the gravitational redshift on the sub-
centimeter scale with 5 spatially multiplexed ensembles of 87Sr. We observe a
gravitational redshift for millimeter to centimeter scale differences in
height, and find that it is consistent with the gravitational redshift
expected from general relativity to within $1\sigma$ total uncertainty. We
demonstrate a gravitational redshift measurement resolution of 1.3 mm. Our
results highlight the use of the spatially multiplexed ensemble techniques for
achieving long coherence times and low differential instabilities without the
need for a state-of-the-art clock laser, and demonstrate its utility for
characterization of spatially varying systematic shifts in optical lattice
clocks on the sub-centimeter scale and at the $10^{-19}$ level. These results
pave the way for relativistic geodesy at the sub-centimeter scale with optical
atomic clocks [9, 11, 10, 12, 8], the realization of entanglement enhanced
clocks [20, 21, 22], and explorations of the interplay between quantum
mechanics and gravity [23, 24].

Figure 1: $|$Experimental system and synchronous differential clock
comparisons. a, A representative camera image of a spatially multiplexed array
with 5 ensembles of 87Sr atoms (indexed 1–5 from top to bottom) trapped in a
vertical 1D optical lattice for differential clock comparisons. The spacing
between neighboring ensembles is 0.25 cm, spanning a total height difference
of 1 cm. Due to the gravitational redshift, clocks at a higher gravitational
potential are predicted to tick faster than clocks at a lower potential. The
grey box shows the orientations of the applied bias magnetic field ($B$), and
the lattice and clock laser polarizations ($\epsilon$). b, A representative
outcome from synchronous Ramsey spectroscopy on the
$\ket{{}^{1}S_{0},m_{F}=+5/2}\leftrightarrow\ket{{}^{3}P_{0},m_{F}=+3/2}$
clock transition with a 10 s free precession time using 5 atomic ensembles,
resulting in 10 pairwise clock comparisons. In each plot, the excitation
fractions of ensemble $j$ are plotted against the excitation fractions of
ensemble $i$ ($P_{j}:P_{i}$, where $j>i$), tracing out an ellipse which is
fitted to in order to extract the frequency difference between that pair. The
frequency difference is dominated by the first-order differential Zeeman
shift, which is rejected when averaging transitions between opposite spin
states (Methods).

Figure 2: $|$Sources and characterization of primary systematic shifts. a, The
differential density shift at $u_{op}=15\leavevmode\nobreak\ E_{\rm rec}$ is
evaluated by varying the atom number difference ($\Delta N$). A linear fit
yields a shift of $-0.7(1)\times 10^{-19}$ per 100 atom number difference. b,
Evaluation of the second-order Zeeman gradient arising from the magnetic field
gradient ($\partial B/\partial z\approx 1.5$ mG$/$cm). A linear fit yields a
gradient of $-95.3(9)\times 10^{-19}/$cm. c, Characterization of the BBR shift
due to thermal gradients across the vacuum chamber. The inset is an
illustration of the science chamber. To evaluate the BBR effect, a thermal
gradient is introduced by varying the temperature difference between top and
bottom viewports by up to $\pm 1$ K. A linear fit yields a BBR sensitivity of
$-4.2(1)\times 10^{-18}/$cm per 1 K difference in our system. d, Correlations
between $\delta u$ and differential lattice light shifts after subtraction of
the residual spatial light shift gradient. e, Evaluation of lattice light
shift gradient at $u_{op}=15\leavevmode\nobreak\ E_{\rm rec}$ after removing
contributions from $\delta u$ shifts. A linear fit yields a gradient of
$-8(1)\times 10^{-20}/E_{\rm rec}/$cm.

Figure 3: $|$Gravitational redshift measurements. a, The measured fractional
frequency gradients after accounting for all systematic corrections over 14
data taking runs (red scatter points), each of duration ranges from 1 to 4
hours. The weighted average (red solid line) yields a measured gradient of
$[-12.4\pm 0.7_{(\rm stat)}\pm 2.4_{(\rm sys)}]\times 10^{-19}/$cm, consistent
with the expected redshift gradient (black solid line). Red (blue) area
represents $\pm 1\sigma$ statistical (total) uncertainty, in which the
statistical uncertainty is scaled up by the square root of the reduced
$\chi^{2}$ statistic, $\chi^{2}_{\rm red}=1.62$. b, Mean differential
frequencies as a function of height difference across all measurements (red
scatter points), analyzed using the same data set as in a. A linear fit (red
solid line) yields a fractional frequency gradient of $(-11.9\pm 2.5)\times
10^{-19}/$cm, again fully consistent with the expected redshift gradient
(black solid line).

Figure 4: $|$Extracting relative height differences using relativistic
geodesy. The relative clock height differences across the array are determined
using the measured gravitational redshifts and the independently measured
local gravitational acceleration $g$. The double arrows represent the
extracted height difference and the associated uncertainty for each clock
pair. The true clock heights are shown on the right (red values), with the
lowest clock (clock 5) defined as being at a height of 0 cm. All units are in
cm unless otherwise specified.

Sources | Gradient | Uncertainty
---|---|---
| ($\times 10^{-19}/$cm) | ($\times 10^{-19}/$cm)
BBRa | -15.7 | 1.5
Lattice light | -11.8 | 1.1
Densityb | – | 1.0
$2^{\rm{nd}}$ order Zeemana | -95.3 | 0.9
Probe Stark | 0 | 0.5
DC Stark | 0 | 0.1
Ellipse fittingb | – | 0.5
Total systematic correction | +122.8 | 2.4
Statistical gradient | -135.2 | 0.7
Corrected gradient | -12.4 | 2.5
Expected redshift gradient | -10.9 |  $<$0.1

Table 1. Fractional frequency gradients and corresponding uncertainties.
Uncertainties are quoted as 1$\sigma$ standard deviations. For each systematic
effect, more discussion can be found in the Methods.
a The BBR shift and second-order Zeeman shift are corrected for each
measurement, here the weighted averaged values across all 14 measurements are
listed.
b The density shift and bias error from ellipse fitting are corrected for each
pairwise clock comparison individually, and do not depend linearly on spatial
separation.

## References

  * [1] Einstein, A. Grundgedanken der allgemeinen Relativitätstheorie und Anwendung dieser Theorie in der Astronomie. Preuss. Akad. der Wissenschaften. Sitzungsberichte 315, 778-786 (1915)
  * [2] Hafele, J. C. & Keating, R. Around-the-world atomic clocks: predicted relativistic time gains. Science 177, 166–168 (1977)
  * [3] Chou, C. W., Hume, D. B., Rosenband, T. & Wineland, D. J. Optical clocks and relativity. Science 329, 1630–1633 (2010)
  * [4] Delva, P. et al. Gravitational redshift test using eccentric Galileo satellites. Phys. Rev. Lett. 121, 231101 (2018)
  * [5] Herrmann, S. et al. Test of the gravitational redshift with Galileo satellites in an eccentric orbit. Phys. Rev. Lett. 121, 231102 (2018)
  * [6] Takamoto, M. et al. Test of general relativity by a pair of transportable optical lattice clocks. Nat. Photon. 14, 411–415 (2020)
  * [7] Bothwell, T. et al. Resolving the gravitational redshift across a millimetre-scale atomic sample. Nature 602, 420–424 (2022)
  * [8] McGrew, W. F. et al. Atomic clock performance enabling geodesy below the centimetre level. Nature 564, 87–90 (2018)
  * [9] Delva, P. & Lodewyck, J. Atomic clocks: new prospects in metrology and geodesy. Acta Futura 7, 67–78 (2013)
  * [10] Lion, G. et al. Determination of a high spatial resolution geopotential model using atomic clock comparisons. J. Geod. 91, 597–611 (2017)
  * [11] Takano, T. et al. Geopotential measurements with synchronously linked optical lattice clocks. Nat. Photon. 10, 662-666 (2016)
  * [12] Grotti, J. et al. Geodesy and metrology with a transportable optical clock. Nat. Phys. 14, 437–441 (2018)
  * [13] Derevianko, A. & Pospelov, M. Hunting for topological dark matter with atomic clocks. Nat. Phys. 10, 933–936 (2014)
  * [14] Safronova, M. S. et al. Search for new physics with atoms and molecules. Rev. Mod. Phys. 90, 025008 (2018)
  * [15] Wsciło, P. et al. New bounds on dark matter coupling from a global network of optical atomic clocks. Sci. Adv. 4, eaau4869 (2018)
  * [16] Kennedy, C. J. et al. Precision metrology meets cosmology: improved constraints on ultralight dark matter from atom-cavity frequency comparisons. Phys. Rev. Lett. 125, 201302 (2020)
  * [17] Delaunay, C., Ozeri, R., Perez, G. & Soreq, Y. Probing atomic Higgs-like forces at the precision frontier. Phys. Rev. D 96, 093001 (2017)
  * [18] Counts, I. & Hur, J. et al. Evidence for nonlinear isotope shift in Yb+ search for new boson. Phys. Rev. Lett. 125, 123002 (2020)
  * [19] Kolkowitz, S. et al. Gravitational wave detection with optical lattice atomic clocks. Phys. Rev. D 94, 124043 (2016)
  * [20] Gil, L. I. R., Mukherjee, R., Bridge, E. M., Jones, M. P. A. & Pohl, T. Spin squeezing in a Rydberg lattice clock. Phys. Rev. Lett. 112, 103601 (2014)
  * [21] Pedrozo-Peñafiel, E. et al. Entanglement on an optical atomic-clock transition. Nature 588, 414–-418 (2020)
  * [22] Van Damme, J., Zheng, X., Saffman, M., Vavilov, M. G. & Kolkowitz, S. Impacts of random filling on spin squeezing via Rydberg dressing in optical clocks. Phys. Rev. A 103, 023106 (2021)
  * [23] Colella, R. and Overhauser, A. W. & Werner, S. A. Observation of gravitationally induced quantum interference. Phys. Rev. Lett. 34, 1472 (1975)
  * [24] Page, Don N. & Geilker, C. D. Indirect evidence for quantum gravity. Phys. Rev. Lett. 47, 979 (1981)
  * [25] Will, C. M. The confrontation between general relativity and experiment. Living Rev. Relativ. 17, 4 (2014)
  * [26] Dyson, F. W., Eddington, A. S. & Davidson, C. IX. A determination of the deflection of light by the sun’s gravitational field, from observations made at the total eclipse of May 29, 1919. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 220, 291–333 (1920)
  * [27] Clemence, G. M. The relativity effect in planetary motions. Rev. Mod. Phys. 19, 361–364 (1947)
  * [28] Pound, R. V. & Rebka, G. A. Gravitational red-shift in nuclear resonance. Phys. Rev. Lett. 3, 439 (1959)
  * [29] Ciufolini, I. & Pavlis, E. C. A confirmation of the general relativistic prediction of the lense–thirring effect. Nature 431, 958–960 (2004)
  * [30] Stairs, I. H. Testing general relativity with pulsar timing. Living Rev. Relativ. 6, 5 (2003)
  * [31] Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration) Observation of gravitational waves from a binary black hole merger. Phys. Rev. Lett. 116, 061102 (2016)
  * [32] Archibald, A. M. et al. Universality of free fall from the orbital motion of a pulsar in a stellar triple system. Nature 559, 73–76 (2018)
  * [33] Ludlow, A. D., Boyd, M. M., Ye, J., Peik, E. & Schmidt, P. O. Optical atomic clocks. Rev. Mod. Phys. 87, 637–701 (2015)
  * [34] Huntemann, N., Sanner, C., Lipphardt, B., Tamm, C. & Peik, E. Single-ion atomic clock with $3\times{10}^{-18}$ systematic uncertainty. Phys. Rev. Lett. 116, 063001 (2016)
  * [35] Bothwell, T. et al. JILA SrI optical lattice clock with uncertainty of $2.0\times 10^{-18}$. Metrologia 56, 065004 (2019)
  * [36] Brewer, S. M. et al. ${}^{27}{\rm Al}^{+}$ quantum-logic clock with a systematic uncertainty below ${10}^{-18}$. Phys. Rev. Lett. 123, 033021 (2019)
  * [37] Oelker, E. et al. Demonstration of $4.8\times 10^{-17}$ stability at 1 s for two independent optical clocks. Nat. Photon. 13, 714–719 (2019)
  * [38] Campbell, S. L. et al. A Fermi-degenerate three-dimensional optical lattice clock. Science 358, 90–94 (2017)
  * [39] Koller, S. B. et al. Transportable optical lattice clock with $7\times 10^{-17}$ uncertainty. Phys. Rev. Lett. 118, 073601 (2017)
  * [40] Kessler, T. et al. A sub-40-mHz-linewidth laser based on a silicon single-crystal optical cavity. Nat. Photon. 6, 687–692 (2012)
  * [41] Zhang, W. et al. Ultrastable silicon cavity in a continuously operating closed-cycle cryostat at 4 K. Phys. Rev. Lett. 119, 243601 (2017)
  * [42] Robinson, J. M. et al Crystalline optical cavity at 4 K with thermal-noise-limited instability and ultralow drift. Optica 6, 240–243 (2019)
  * [43] Clements, E. R. et al. Lifetime-limited interrogation of two independent ${}^{27}{\mathrm{Al}}^{+}$ clocks using correlation spectroscopy. Phys. Rev. Lett. 125, 243602 (2020)
  * [44] Schioppo, M. et al. Ultrastable optical clock with two cold-atom ensembles. Nat. Photon. 11, 48–52 (2017)
  * [45] Zheng, X. et al. Differential clock comparisons with a multiplexed optical lattice clock. Nature 602, 425–430 (2022)
  * [46] Kim, M. E. et al. Optical coherence between atomic species at the second scale: improved clock comparisons via differential spectroscopy. Preprint at https://arxiv.org/abs/2109.09540 (2021)
  * [47] Brown, R. C. et al. Hyperpolarizability and operational magic wavelength in an optical lattice clock. Phys. Rev. Lett. 119, 253001 (2017)
  * [48] Ushijima, I., Takamoto, M. & Katori, H. Operational magic intensity for Sr optical lattice clocks. Phys. Rev. Lett. 121, 263202 (2018)
  * [49] Martin, M. J. et al. A quantum many-body spin system in an optical lattice clock. Science 341, 632–636 (2013)
  * [50] Zhang, X. et al. Spectroscopic observation of SU($N$)-symmetric interactions in Sr orbital magnetism. Science 345, 1467–1473 (2014)
  * [51] Aeppli, A. et al. Hamiltonian engineering of spin-orbit coupled fermions in a Wannier-Stark optical lattice clock. Preprint at https://arxiv.org/abs/2201.05909 (2022)

##  Methods

###  Sample preparation and experimental procedure

The experimental sequence starts with laser cooling the atoms down to
$1\leavevmode\nobreak\ \mu$K temperature with standard two-stage magneto-
optical trapping (MOT). Using the multiplexed ensemble loading technique with
a movable one-dimensional (1D) optical lattice described and demonstrated in
Ref. [45], 5 ensembles of ultra-cold, spin-mixed 87Sr atoms are loaded into a
vertical optical lattice with a depth of $130\leavevmode\nobreak\
E_{\textrm{rec}}$, where $E_{\textrm{rec}}/h\approx 3.5$ kHz is the recoil
energy of a lattice photon and $h$ is the Planck constant, with an equal
spacing between ensembles of 0.25 cm over a total extent of 1 cm vertically.
This is followed by hyperfine spin polarization into either stretched state
($\ket{{}^{1}S_{0},m_{F}=\pm 9/2}$) and in-lattice cooling. The lattice is
then adiabatically ramped down to the operational trap depth
($u_{op}=15\leavevmode\nobreak\ E_{\textrm{rec}}$), at which a series of $\pi$
pulses addressing the ${}^{1}S_{0}\leftrightarrow{}^{3}P_{0}$
($g\leftrightarrow e$) transitions prepare the atoms into $\ket{e,m_{F}=\pm
3/2}$. Ramsey spectroscopy is performed by interrogating the $\ket{g,m_{F}=\pm
5/2}\leftrightarrow\ket{e,m_{F}=\pm 3/2}$ clock transitions. Following the
second Ramsey $\pi/2$ pulse, the lattice is adiabatically ramped back up to
130 $E_{\textrm{rec}}$ for read-out. The populations in the ground and excited
clock states of all 5 ensembles are read-out in parallel with imaging pulses
along the lattice axis, with scattered photons collected on a camera (Andor,
iXon-888). The excitation fraction is extracted through
$P=(N_{e}-N_{bg})/(N_{g}+N_{e}-2N_{bg})$, where $N_{g},N_{e}$ and $N_{bg}$ are
the ground state population, excited clock state population, and background
counts without atoms, respectively.

As described in Ref. [45], for multiple ensemble preparation, we chirp the
frequency of the retro-reflected lattice beam during the single-frequency
stage of the narrow-line ${}^{1}S_{0}\leftrightarrow{}^{3}P_{1}$ second-stage
MOT. We perform 4 lattice movements of 0.25 cm each in order to load 5
spatially separated atom ensembles, and a final lattice movement of 0.5 cm in
the opposite direction, which positions the ensemble array symmetrically
around the lattice beam waist. The entire duration of the moving lattice
portion of the loading sequence is typically less than 100 ms, with 80% of
total atom number transferred efficiently from the narrow-line MOT to the 5
ensembles. In-lattice axial (sideband) and radial (Doppler) cooling are
subsequently applied to lower the atom temperatures after lattice
acceleration.

For clock interrogation, we probe the $\ket{g,m_{F}=\pm
5/2}\leftrightarrow\ket{e,m_{F}=\pm 3/2}$ transitions with a shared clock
laser along the lattice axis. Synchronous Ramsey spectroscopy is performed to
reject common-mode local oscillator noise. The typical Ramsey interrogation
time ($T_{R}$) is roughly 10 s, with a dead time ($T_{d}$) of 2 s between
interrogations for sample preparation and read-out, yielding a measurement
duty cycle of 83%. Simultaneously probing 5 ensembles results in 10 pairwise
clock comparisons for a single nuclear spin state. Combined with atom numbers
($N$) of $\sim 2000$ per ensemble and contrast ($C$) of above 80 %, the
typical differential instability for each pairwise comparison is below
$1\times 10^{-17}/\sqrt{\tau}$, where $\tau$ is the averaging time, consistent
with the quantum projection noise (QPN) limit

$\sigma_{\rm QPN}(\tau)=\frac{\sqrt{2}}{2\pi
fCT_{R}}\sqrt{\frac{T_{R}+T_{d}}{N\tau}},$ (3)

where $f$ is the clock frequency, $\tau$ is the averaging time, and the factor
of $\sqrt{2}$ assumes equal contribution from each clock.

###  Atomic temperatures

Representative data from clock sideband thermometry on the
$\ket{g}\leftrightarrow\ket{e}$ transition at $130\leavevmode\nobreak\
E_{\textrm{rec}}$ is shown in Extended Data Fig. 1a. The reduced height of the
red sideband ($n_{z}\rightarrow n_{z}-1$, where $n_{z}$ is the axial
vibrational quantum number) indicates that $>99\%$ of the atoms are populated
in the lattice ground band ($n_{z}=0$). A separate clock probe perpendicular
to the lattice axis taken at 15 $E_{\textrm{rec}}$ is shown in Extended Data
Fig. 1b, in which the Doppler broadened profile indicates a radial atomic
temperature below $200$ nK.

###  Atomic coherence times

We use a synchronized Ramsey sequence to measure the atomic coherence. We
prepare two atomic ensembles in the $\ket{e,m_{F}=+3/2}$ state and apply a
$\pi$/2 pulse on the $\ket{e,m_{F}=+3/2}\leftrightarrow\ket{g,m_{F}=+5/2}$
transition. After waiting for a varying dark time ranging from 100 ms to 40 s,
we apply a second $\pi$/2 pulse with a random phase relative to the first. We
then measure the excitation fractions of each ensemble and plot the averaged
contrast as a function of dark times. We observe an exponential decay of the
contrast with a time constant of $32.1(9)$ s at $15\leavevmode\nobreak\ E_{\rm
rec}$ (Extended Data Fig. 2), 300 times longer than the measured atom-laser
coherence time ($\sim 100$ ms). We note that while the atomic coherence time
is about 30 s, the typical clock interrogation time in our measurements is
around 10 s, chosen to minimize the QPN because lost of atoms and reduced
contrast at longer dark times [45].

###  Ellipse phase extraction

We perform synchronous Ramsey spectroscopy with up to 5 ensembles (indexing
$1,2,\dots,5$ from top to bottom). This results in 10 pairs of clock
comparisons performed simultaneously when probing transition from either
nuclear spin states $\ket{g,m_{F}=\pm 5/2}$. For each pair of clock
comparison, we plot the excitation fractions $P_{i}(P_{j})$ of ensemble $i(j)$
on the x(y)-axis (note that we choose the convention $j>i$). The excitation
fractions are given by

$\begin{split}P_{i}&=\frac{1}{2}\bigg{[}1+C_{i}{\rm
cos}(\theta_{L})\bigg{]},\\\ P_{j}&=\frac{1}{2}\bigg{[}1+C_{j}{\rm
cos}(\theta_{L}+\phi_{d})\bigg{]},\end{split}$ (4)

where $C_{i(j)}$ is the contrast of ensemble $i(j)$, $\theta_{L}$ is the
common-mode laser phase, and $\phi_{d}$ is the differential phase which yields
the differential frequency ($\delta f_{ij}=f_{j}-f_{i}$) between ensemble pair
$(i,j)$ through $\phi_{d}=2\pi\delta f_{ij}T_{R}$ for a given known Ramsey
free evolution time $T_{R}$.

Since we are operating at Ramsey dark times well beyond the laser coherence
times, $\theta_{L}$ is random and uniformly distributed from $0$ to $2\pi$.
The data randomly samples from points lying on an ellipse (with slight
deviations from the ellipse due to QPN). We then fit to this ellipse using a
least-squares approach [1]. To extract the differential phase [2], we rewrite
the data $\\{P_{i},P_{j}\\}$ (denoted as $\\{x,y\\}$ below) in the form of a
generalized conic section

$a_{1}x^{2}+a_{2}xy+a_{3}y^{2}+a_{4}x+a_{5}y+a_{6}=0,$ (5)

which describes an ellipse when $a_{2}^{2}-4a_{1}a_{3}<0$. We rewrite Eq. 4 as

$\begin{split}x^{\prime}&=\frac{2x-1}{C_{x}}={\rm cos}(\theta_{L}),\\\
y^{\prime}&=\frac{2y-1}{C_{y}}={\rm cos}(\theta_{L}+\phi_{d}).\end{split}$ (6)

Through cancelling out $\theta_{L}$, we have

$\begin{split}&\frac{4}{C_{x}^{2}}x^{2}-\frac{8{\rm
cos}\phi_{d}}{C_{x}C_{y}}+\frac{4}{C_{y}^{2}}y^{2}+\bigg{(}\frac{4}{C_{x}C_{y}}-\frac{4}{C_{x}^{2}}\bigg{)}x+\bigg{(}\frac{4}{C_{x}C_{y}}-\frac{4}{C_{y}^{2}}\bigg{)}y\\\
&+\bigg{(}\frac{1}{C_{x}^{2}}-\frac{2}{C_{x}C_{y}}+\frac{1}{C_{y}^{2}}-{\rm
sin}^{2}\phi_{d}\bigg{)}=0,\end{split}$ (7)

which can be matched up with the coefficients $a_{i}$ from Eq. 5. The
differential phase $\phi_{d}$ is then extracted using:

$\phi_{d}={\rm cos}^{-1}\bigg{(}\frac{-a_{2}}{2\sqrt{a_{1}a_{3}}}\bigg{)}.$
(8)

The associated Allan deviation is extracted via jackknifing technique [3], and
is then fitted to a white frequency noise model with $1/\sqrt{\tau}$ scaling.
Extrapolating the fit to the full averaging time yields the statistical
uncertainty of the differential frequency.

In our measurements, we probe with an interleaved sequence between clock
transitions with either nuclear spin state, $\ket{g,m_{F}=\pm
5/2}\leftrightarrow\ket{e,m_{F}=\pm 3/2}$. This results in 10 ellipses for
transition with a single nuclear spin state, and thus 20 ellipses per
measurement. A representative plot of the
$\ket{g,m_{F}=+5/2}\leftrightarrow\ket{e,m_{F}=+3/2}$ transition is shown in
Fig. 1b. We note that the differential phase for each ellipse is dominated by
the differential first-order Zeeman shift (see section “Systematic
evaluation”), which is on the order of $\pm 8\times 10^{-17}/$cm and is
rejected by averaging transitions with opposite spin states.

###  Data blinding protocol

To eliminate possible bias of our data taking and systematic analysis towards
an expected outcome, we employ a data blinding protocol. Our data software
adds a large constant offset gradient to our measurements, including the data
taken for systematic evaluations and data runs taken under normal operating
conditions. The blinded offset gradient is pseudo-randomly drawn from a
uniform distribution spanning over $\pm 5\times 10^{-18}/$cm, 10 times the
size of the expected redshift gradient. The blinded offset gradient is scaled
by the height difference between each ensemble pair, and is then automatically
added to the results from our data analysis code for ellipse phase extraction.

The blinded offset gradient was only unblinded after finalizing the
corrections for all systematic effects, determining the measured value with
blinded offset gradient taken under normal operating conditions, and
finalizing the associated statistical and systematic uncertainties. No
additional data was taken and no changes were made to the analysis, the error
budget, the measured value, or the uncertainties after unblinding.

###  Systematic evaluation

####  Zeeman shifts

Coupling of the clock states to external magnetic fields give rise to Zeeman
shifts of the clock transition frequency. For a single atomic ensemble with
magnetic field amplitude $B$, the Zeeman shifts when probing the
$\ket{g,m_{F}=\pm 5/2}\leftrightarrow\ket{e,m_{F}=\pm 3/2}$ clock transitions
can be expressed as

$\nu_{\rm ZS}(B,\pm)=\pm\mu_{L}B+\mu_{Q}B^{2},$ (9)

where $\mu_{L}$ and $\mu_{Q}$ are the first and second-order Zeeman shift
coefficients, respectively. For differential clock comparison between ensemble
pair ($i,j$), where $j>i$, with differing magnetic field strengths $B_{i}$ and
$B_{j}$, the differential Zeeman shifts, $\Delta\nu_{\rm
ZS}(B_{i},B_{j},\pm)=\nu_{\rm ZS}(B_{j},\pm)-\nu_{\rm ZS}(B_{i},\pm)$, are
given by

$\begin{split}\Delta\nu_{\rm
ZS}(B_{i},B_{j},\pm)&=\pm\mu_{L}(B_{j}-B_{i})+\mu_{Q}(B_{j}^{2}-B_{i}^{2}),\\\
&=\pm\mu_{L}\delta B+2\mu_{Q}B\delta B,\end{split}$ (10)

where $B=(B_{i}+B_{j})/2$ is the mean magnetic field amplitude, and $\delta
B=B_{j}-B_{i}$ is the magnetic field difference. The first-order differential
Zeeman shift (on the order of $\pm 8\times 10^{-17}/$cm) is rejected by
averaging opposite spin state transitions, while their splitting yields
$\delta B$. In the limit of $|\delta B|\ll B$, the second-order differential
Zeeman shift can be approximated as

$\begin{split}\delta B&=\Delta\nu_{\rm split}/(2\mu_{L}),\\\ \Delta\nu_{\rm
ZS,2}&=\xi\Delta\nu_{\rm split}B,\end{split}$ (11)

where $\xi=\mu_{Q}/\mu_{L}$, and $\Delta\nu_{\rm split}$ is the splitting
between the transitions with opposite spin states after subtraction of
residual vector AC Stark shift.

The second-order Zeeman shift correction is applied for each measurement run,
with uncertainty primarily limited by $\xi$, which is found to be $0.0105(1)$
G-1 by varying $B$ and $\delta B$. This is in good agreement with the
theoretical value of $0.0104(2)$, calculated using
$\mu_{L}=-22.4\leavevmode\nobreak\ $Hz$/$G and $\mu_{Q}=-0.233(5)$ Hz$/$G2
[4]. The bias magnetic field $B$ is calibrated each measurement day using the
magnetically sensitive transitions $\ket{{}^{1}S_{0},m_{F}=\pm
9/2}\rightarrow\ket{{}^{3}P_{1},m_{F}=\pm 11/2}$, and has a typical value of
$\sim 5.5$ G with fractional uncertainty below $10^{-4}$. The field difference
$\delta B$ ($\sim 1.5$ mG$/$cm) is extracted from the splitting between the
clock transitions with opposite nuclear spin states with fractional
uncertainty below $10^{-3}$ after subtraction of the vector AC Stark shift.
Through weighted averaging of the data taken at normal operating conditions,
we find the mean second-order Zeeman gradient to be $-95.3(9)\times
10^{-19}/$cm (Fig. 2b).

####  Density shift

Due to the Pauli-exclusion principle, $s$-wave interactions are forbidden for
identical Fermionic atoms within a single lattice site, while $p$-wave
collisions are allowed, leading to a clock frequency shift that scales
linearly with atomic density [49, 50]. In our system, the differential density
shift is evaluated by varying the atom numbers loaded into each ensemble. We
find a linear slope of $-0.7(1)\times 10^{-19}$ per 100 atom number difference
at our operational lattice trap depth $u_{op}=15\leavevmode\nobreak\ E_{\rm
rec}$ (Fig. 2a). For a typical run of the experiment, each atomic ensemble has
about 2000 atoms, corresponding to an overall density shift of about
$-14(2)\times 10^{-19}$. The differential density shift is suppressed by a
factor of 10 when the atom number difference is bounded below 200 through
optimization of ensemble loading times and conditions. In addition, we measure
the atom number in each ensemble in every shot of the experiment, and the
average differential density shift is calculated and corrected individually
for each pairwise clock comparison measurement run. Uncertainty of density
shift correction arises from shot-to-shot atom number fluctuations, for which
we estimate an upper bound limit of $1\times 10^{-19}$.

We note that a recent study found that operation in a gravity-tilted shallow
lattice allows for cancellation of density shifts at a “magic” lattice depth
near $12\leavevmode\nobreak\ E_{\rm rec}$, where the partially delocalized
Wannier-Stark states enable tunability of on-site $p$-wave versus
neighbouring-site $s$-wave atomic interactions [51]. In this work we did not
observe such a cancellation effect at shallower lattice depths, likely due to
the difference in dynamics between long interrogation time Rabi spectroscopy
as was done in Ref. [51] versus our use of Ramsey spectroscopy with a 50:50
superposition [5]. This offers the prospect of further reducing uncertainty
from differential density shifts in future works.

####  Black body radiation shift

In our experiment, the optical lattice is orientated nearly vertically, with a
tilt of about $5^{\circ}$, and is centered with respect to the science chamber
to the best of abilities. The tilt is determined by measuring the Wannier-
Stark ladder resonances (refer to “Wannier-Stark ladder and the expected
gravitational redshift” for details). The recessed high-emissivity Fused
Silica viewports mounted on the top and bottom sides of the science chamber
are the closest surfaces to the atoms, and are primarily responsible for the
black body radiation (BBR) gradient along the lattice axis. To study the BBR
effect [6, 7], we heat up either the top or bottom stainless steel flange of
the science chamber, which results in a temperature difference of up to $\pm
1$ K between the top and bottom viewports (see inset of Fig. 2c). By measuring
the frequency differences of the 10 ensemble pairs simultaneously, we find
that the resulting frequency shifts scale linearly with the temperature
difference, and the corresponding slopes scale roughly linearly with the
height differences (Extended Data Fig. 3). Through a linear fit to the
extracted slopes as a function of height differences, we find the BBR
sensitivity in our system to be $-4.2(1)\times 10^{-18}/$cm per 1 K difference
(Fig. 2c).

To monitor the temperature, we use commercially available 10 k$\Omega$
negative temperature coefficient thermistors (Amphenol Thermometrics,
MC65F103A), rated for an interchangeability of 50 mK by the manufacturer.
Relative accuracy of the temperature sensors is calibrated in an ice water
bath, and we find the temperature differences between each sensor consistent
within 25 mK. The temperature sensors are mounted on the stainless steel
Conflat flanges of the top and bottom viewports, as well as the side flanges
of the science chamber for monitoring the radial temperature inhomogeneity.
Under normal operating conditions, we find the temperature difference between
the top and bottom sensors lies within a range of $350\pm 100$ mK and find the
radial temperature inhomogeneity to be less than 75 mK.

We note that the MOT coils are mounted in the recessed window of the science
chamber, which could potentially introduce temporal thermal drifts during the
sample preparation stage. To suppress this, we opt for compact coils with
efficient water cooling through hollow wires, keep the
${}^{1}S_{0}\leftrightarrow{}^{1}P_{1}$ first-stage MOT loading time ($<0.5$
s) short compared to the experimental cycle time ($12$ s), and ensure that the
duty-cycles remain consistent throughout the measurements.

While the strontium oven atom source (AOSense) is another potential source of
BBR shift, the atom ensembles have no direct line of sight to the BBR photons
from the oven thanks to the use of a 2D-MOT for atomic beam deflection. We
operate the oven at the lowest possible temperature
($T=360\leavevmode\nobreak\ ^{\circ}$C) under normal conditions, and find no
statistically significant change in the frequency differences across the
ensemble array at the $1\times 10^{-19}/$cm level when intentionally
increasing the oven temperature to 460 ∘C, a lever arm of $1.8$ due to $T^{4}$
scaling bounds the systematic gradient from the oven to below $5.5\times
10^{-20}/$cm.

Overall, for measurements taken under normal operations, we find the averaged
temperature difference between top and bottom sensors to be $370(35)$ mK,
resulting in a BBR gradient of $-15.7(1.5)\times 10^{-19}/$cm.

####  Lattice light shift

The lattice light shift for a single atomic ensemble is given by [48]:

$\begin{split}\nu_{LS}(u,\delta_{L},n_{z})\approx&\bigg{(}\frac{\partial\tilde{\alpha}^{E1}}{\partial\nu}\delta_{L}-\tilde{\alpha}^{qm}\bigg{)}\bigg{(}n_{z}+\frac{1}{2}\bigg{)}u^{1/2}\\\
&-\bigg{[}\frac{\partial\tilde{\alpha}^{E1}}{\partial\nu}\delta_{L}+\frac{3}{2}\tilde{\beta}\bigg{(}n_{z}^{2}+n_{z}+\frac{1}{2}\bigg{)}\bigg{]}u\\\
&+2\tilde{\beta}\bigg{(}n_{z}+\frac{1}{2}\bigg{)}u^{3/2}-\tilde{\beta}u^{2},\end{split}$
(12)

in which $\tilde{\alpha}^{E1}$, $\tilde{\alpha}^{qm}$, and $\tilde{\beta}$ are
the differential $E$1, $E$2-$M$1 polarizabilities, and hyperpolarizability on
the clock transition, respectively. $u$ is the lattice trap depth in units of
$E_{\textrm{rec}}$, $n_{z}$ is the axial vibrational quanta, and $\delta_{L}$
is lattice detuning from the effective magic wavelength, where the scalar and
tensor shifts cancel. The vector shift is rejected by averaging the
transitions with opposite nuclear spin states. Through clock sideband
thermometry both axially and radially we find $n_{z}<0.01$ and radial
temperature $T_{r}<200$ nK at $u_{op}=15\leavevmode\nobreak\ E_{\rm rec}$,
where thermal averaging of the effective trap depth can be neglected. By
operating at lattice depths below $35\leavevmode\nobreak\ E_{\rm{rec}}$, the
hyperpolarizability terms are also negligible.

To model the differential lattice light shift in our system, we introduce a
dimensionless parameter $\delta u$, which characterizes the relative lattice
trap depth difference between ensemble pair $(i,j)$, where $j>i$

$\begin{split}u_{i}&=u,\\\ u_{j}&=u(1+\delta u).\end{split}$ (13)

For $\delta u\ll 1$, the differential light shift ($\nu_{LS,j}-\nu_{LS,i}$)
can be approximated as

$\Delta\nu_{LS}(\delta_{L},\delta u,u)=\delta
u\bigg{[}\bigg{(}\frac{\partial\tilde{\alpha}^{E1}}{\partial\nu}\delta_{L}-\tilde{\alpha}^{qm}\bigg{)}\frac{1}{4}\sqrt{u}-\frac{\partial\tilde{\alpha}^{E1}}{\partial\nu}\delta_{L}u\bigg{]},$
(14)

which scales linearly with $\delta u$.

By modulating the lattice detuning between $\delta_{L}+100$ MHz and
$\delta_{L}-100$ MHz at the operational depth $u_{op}=15\leavevmode\nobreak\
E_{\rm rec}$, we have

$\begin{split}&\Delta\nu_{LS}(u_{op},\delta_{L}+100\leavevmode\nobreak\ {\rm
MHz})-\Delta\nu_{LS}(u_{op},\delta_{L}-100\leavevmode\nobreak\ {\rm MHz})\\\
=&\leavevmode\nobreak\ \delta
u\frac{\partial\tilde{\alpha}^{E1}}{\partial\nu}\bigg{(}\frac{1}{4}\sqrt{u_{op}}-u_{op}\bigg{)}\times
200\leavevmode\nobreak\ {\rm MHz},\end{split}$ (15)

allowing extraction of $\delta u$, which are found to be below 5 % for all the
ensemble pairs, and are symmetric around 0 (Extended Data Fig. 4a), reflecting
the Gaussian nature of the lattice beam profile and that the array is centered
about the focus. Through mapping out the frequency shifts for 10 ensemble
pairs when modulating the lattice intensity, we find the differential light
shifts scale with the spatial separations between ensemble pairs but do not
scale with $\delta u$, resulting in a residual spatial light shift gradient of
$-8(1)\times 10^{-20}/E_{\rm rec}/$cm (Extended Data Fig. 4b).

Although the origin of the spatial light shift gradient is not definitively
known, we also observe a differential lattice vector Stark gradient of
$-2.5(2)\times 10^{-18}/E_{\rm rec}/$cm by measuring the splitting between
transitions with opposite spin states (Extended Data Fig. 4c), which also does
not scale with $\delta u$ and instead scales with trap depth and spatial
separation. For a single ensemble, the vector Stark shift is proportional to
$\xi_{\rm lat}(\hat{e}_{k}\cdot\hat{e}_{B})$, where $\xi_{\rm lat}$ is the
ellipticity of the lattice light, $\hat{e}_{k}$ is the lattice wave vector
along $\hat{z}$ direction, and $\hat{e}_{B}$ is the magnetic field vector
primarily along $\hat{x}$. We find the vector Stark gradient arises from the
spatially varying coefficient, $\hat{e}_{k}\cdot\hat{e}_{B}$, verified through
the change of vector Stark gradient by applying an additional magnetic field
gradient $\partial B/\partial z$ of up to $\pm 10$ mG$/$cm, which effectively
changes $(\hat{e}_{k}\cdot\hat{e}_{B,j}-\hat{e}_{k}\cdot\hat{e}_{B,i})$. Based
on this observation, we hypothesize that the residual spatial light shift
gradient is due to a differential tensor Stark shift, in which the spatially
varying magnetic field vectors are coupled to the (nearly) linear lattice
polarization ($\hat{e}_{\rm lat}$, along $\hat{x}$) through $|\hat{e}_{\rm
lat}\cdot\hat{e}_{B}|^{2}$. We note that the differential vector Stark shift
results in a gradient of $-37(3)\times 10^{-18}/$cm at
$u_{op}=15\leavevmode\nobreak\ E_{\rm rec}$ in the splitting between
transitions with opposite spin states, equivalent to a fictitious magnetic
field gradient of $\sim 600\leavevmode\nobreak\ \mu$G$/$cm and is accounted
for during the second-order Zeeman shift corrections.

To independently evaluate the light shift gradient, we need to account for the
$\delta u$ dependent shifts, which require the knowledge of $\delta u$ and
$\delta_{L}$. The latter is done by first subtracting the residual spatial
gradient from the measured differential light shifts. We then find that the
remaining shifts are correlated with the extracted $\delta u$ (Fig. 2d) as
expected from Eq. 14. From this we find our operational lattice detuning
$\delta_{L}$ to be $16.9(1.5)$ MHz. With the extracted $\delta_{L}$ and
$\delta u$ calibrated from each measurement day, we are able to monitor and
account for the residual light shift gradient. A representative plot is shown
in Fig. 2e.

Our operational lattice frequency is measured to be $368,554,780(30)$ MHz,
limited by the accuracy of the wave-meter (HighFinesse, WS7). However, the
lattice frequency is stabilized to a ultra-low-expansion cavity with typical
drifts of less than $-20$ kHz per day, inferred from the narrow-linewidth
${}^{1}S_{0}\leftrightarrow{}^{3}P_{1}$ MOT, providing sufficient long-term
stability throughout the measurements. The lattice intensity is actively
stabilized and controlled via feedback to the acoustic-optical modulator
before the beam delivering optical fiber (NKT Photonics, LMA-PM-15). Both the
incoming and retro-reflected lattice alignments are monitored on several
cameras and photo-diodes, ensuring the daily-calibrated $\delta u$ remained
symmetric around 0 throughout the data taking campaign. Overall, the lattice
light shift gradient in our system is evaluated to be $-11.8(1.1)\times
10^{-19}/$cm.

####  Probe Stark shift

The probe AC Stark shift arises from the clock light itself, and is suppressed
when probing with a shared clock light, with uncertainties primarily arising
from inhomogeneity and misalignment of the clock beam with respect to the
lattice. The clock light beam waist is $\sim 1$ mm, a factor of 10 greater
than that of the lattice, ensuring homogeneity across the atomic ensembles
both axially and radially. The clock light is carefully aligned to the lattice
light, and is monitored using several cameras by picking-off the beam before
and after the science chamber over a distance of $\sim 2.5$ m. The alignment
is further verified by ensuring that the time periods of the clock transition
Rabi oscillations of all 5 ensembles agree within $1\%$ fitting uncertainty
before and after each measurement run.

Evaluation of the probe Stark effect is performed by interleaving between the
operational clock intensity of $I_{op}$, corresponding to a $\pi$ pulse
duration of 6 ms, and the case of $4I_{op}$. This interleaved measurement is
repeated by further misaligning the clock beam by up to $\pm
500\leavevmode\nobreak\ \mu$rad (estimated via camera images), at which the
$\pi$ pulse durations increase by about a factor of 2 (Extended Data Fig. 5).
Through linear fitting, we find a frequency gradient of $2.5(8)\times
10^{-19}/$cm across the span of $\pm 500\leavevmode\nobreak\ \mu$rad
misalignments. Under normal operations, the misalignment is monitored and
bounded within a $\pm 50\leavevmode\nobreak\ \mu$rad range which accounts for
possible drifts during the experiments. This bounds the uncertainty from probe
Stark shift to below $0.5\times 10^{-19}/$cm.

####  DC Stark shift

Frequency shifts arising from electric fields can perturb the clock transition
through DC Stark effect of the form $kE^{2}$, where $E={\rm|\hat{E}|}$ is the
static electric field and $k\equiv-(\alpha_{e}-\alpha_{g})/(2h)$ is the
coefficient specific to the clock transition, with $\alpha_{g,e}$ being the
static polarizabilities of the ground and excited clock states. To evaluate
the DC Stark effect, a pair of quadrant electrodes are mounted along the
lattice axis outside the top and bottom science chamber viewports, with a
total separation of roughly 30 cm. We then probe the clock transition with
opposite voltages applied to the electrode pair and interleave between
($+V,-V$) and ($-V,+V$) configurations. We observe differential shifts below
the $1\times 10^{-19}/$cm level when applying voltages of up to $\pm 100$ V to
the electrodes, indicating that the background electric field gradient is
small.

To quantify our uncertainty in the differential shifts due to the background
electric field gradient, we follow the approach laid out in Ref. [8]. In the
presence of a background field $E_{bg}$, which would most likely arise due to
charge accumulation on the nearby top and bottom viewports, the DC Stark shift
for a single atom ensemble is given by $k(E_{bg}+cV)^{2}$, where $c$ is the
atomic coefficient. We rewrite the above form as [8]

$\nu_{dc}(V)=\nu_{bg}+aV+bV^{2},$ (16)

where $\nu_{bg}$ is the background stray field shift, and coefficients $a$ and
$b$ are experimentally accessible parameters by modulating $V$. We note that
$aV$ represents the coupling between the background field and the applied
field. When comparing ensemble $i$ and ensemble $j$ ($j>i$), the differential
DC Stark shift ($\Delta\nu_{dc}=\nu_{j}-\nu_{i}$) becomes

$\Delta\nu_{dc}(V)=\Delta\nu_{bg}+AV+BV^{2},$ (17)

where $A=a_{j}-a_{i}$, $B=b_{j}-b_{i}$, and
$\Delta\nu_{bg}=\nu_{bg,j}-\nu_{bg,i}$ is the differential background shift
due to charges on the viewports.

Due to the finite spatial extent ($s$) of the atom ensembles and inhomogeneity
of the applied field, there is no $V$ such that the applied field identically
cancels $\Delta\nu_{bg}$. We consider the extremum value of
$\Delta\nu_{dc}(V)$, denoted as $\Delta\nu^{*}$. We then have the difference
between the differential background shift and the extremum shift,
$\Delta\nu_{\rm diff}\equiv\Delta\nu_{bg}-\Delta\nu^{*}$ (refer to Fig. 1 in
Ref. [8] for details). From Eq. 17, we find $\Delta\nu_{\rm diff}=A^{2}/(4B)$.
We note that $\Delta\nu^{*}$ plays the role of a frequency correction for the
field gradients, and can be written as $\Delta\nu^{*}=\eta\Delta\nu_{\rm
diff}$, due to the fact that $\Delta\nu^{*}$ and $\Delta\nu_{\rm diff}$ scale
similarly with the stray field. The differential background shift is then
given by

$\Delta\nu_{bg}=(1+\eta)\Delta\nu_{\rm diff}=(1+\eta)A^{2}/(4B),$ (18)

in which $\eta$ is given by $\eta=\zeta^{2}s^{2}/R^{2}$ to leading order of
$s$, where $\zeta\equiv(q_{1}+q_{2})/(q_{1}-q_{2})$ quantifies the charge
symmetry between the viewports, and $R$ is the effect length.

By applying voltages of up to $\pm 100$ V to the electrodes, we find
$A^{2}/(4B)=0.6(3)\times 10^{-20}$ for the ensemble pair of 1 cm separation.
In our system, we obtain $\eta\approx 0.02$ for $s=\leavevmode\nobreak\
500\leavevmode\nobreak\ \mu$m, $R=40$ mm, and assuming $25\%$ more charge on
one viewport than the other. This bounds the background DC Stark shift
gradient to below $0.1\times 10^{-19}/$cm.

####  Ellipse fitting

In the presence of QPN and numerical constraints which ensure an ellipse-
specific solution, the least squares ellipse fitting approach is biased at
phases close to 0 or $\pi$ where the ellipse collapses into a straight line
[1]. The differential phases across the ensemble array are dominated by the
differential Zeeman shifts, such that through precise control of the magnetic
field and gradient, we can operate in the regime where all the differential
phases lie within the $(\pi/6,5\pi/6)$ range. For typical experiments with a
10 s Ramsey dark time and 2000 atoms per ensemble, the ellipse fitting bias
error is bound to below $2\times 10^{-19}$ (Extended Data Fig. 6). We note
that during initial data processing, the bias error is corrected using Monte-
Carlo simulations with the experiment parameters such as atom number and
contrast as inputs. While the bias error does not scale with spatial
separations, we estimate an upper bound limit of $0.5\times 10^{-19}$ due to
uncertainty in the input parameters used for Monte-Carlo simulation during
bias error corrections.

####  Wannier-Stark ladder, imaging, and the expected gravitational redshift

The local gravitational acceleration in our laboratory was measured by the
group of Prof. Tikoff from the Department of Geoscience at the University of
Wisconsin-Madison using a LaCoste-Romberg gravimeter, and was then cross-
validated with known values from several other survey points in Wisconsin. The
local gravitational acceleration in our laboratory was measured to be
$g=-9.803$ m$/$s2, rounded to the 4th digit.

The lattice tilt ($\theta_{\rm tilt}$) with respect to gravity is
independently measured using the splitting of the $\pm 1^{\rm st}$ order
Wannier-Stark sidebands [9] to the clock transition at $5\leavevmode\nobreak\
E_{\rm rec}$ lattice depth (Extended Data Fig. 7a), which is given by

$\Delta\nu_{\rm WS,\pm 1}=mg\lambda_{L}{\rm cos}\theta_{\rm tilt},$ (19)

where $m$ is the mass of 87Sr, and $\lambda_{L}=813.4$ nm is the lattice
wavelength. A representative plot of the measured Wannier-Stark sidebands is
shown in Extended Data Fig. 7b, and we find the weighted average splitting
between $\pm 1^{\rm st}$ order sidebands to be 1729(2) Hz, corresponding to a
tilt with respect to gravity of $5.0(1)^{\circ}$.

The overall spatial extent of the ensemble array is $1.00(1)$ cm, calculated
based on the frequency chirp profile of the moving optical lattice loading
sequence, which is controlled using a direct digital synthesizer (DDS, Moglabs
XRF). This yields an effective height difference of $\Delta h=0.99(1)$ cm,
consistent with the extraction of the height difference from the camera
images. The effective camera pixel size in our imaging system is calibrated to
be $34(1)\leavevmode\nobreak\ \mu$m per pixel using the standard time-of-
flight imaging.

Fluctuations of background magnetic fields modifying the narrow-line MOT
operation could result in position drift in the atom ensemble. This is
rejected by identifying the center-of-mass of each ensemble with a 2D Gaussian
fit, followed by a selection of region-of-interest (ROI) with a 15 pixel by 15
pixel region centered on the fit. We note that each ensemble has a finite
vertical spatial extent of roughly $500(10)\leavevmode\nobreak\ \mu$m,
corresponding to an expected gradient of $-5.5(1)\times 10^{-20}$ across the
ensemble due to the gravitational redshift, which is averaged out by spatial
averaging over the ROI in order to extract the redshifts between the center-
of-mass positions of each ensemble. The center-of-mass position uncertainty of
each ensemble is limited by the DDS timing accuracy ($100$ ns), corresponding
a position error of $<0.1\leavevmode\nobreak\ \mu$m.

The expected gravitational redshift from theory of general relativity at a 1
cm height difference is then given by

$\frac{g\Delta h}{c^{2}}=-10.9\times 10^{-19},$ (20)

where $c$ is the speed of light. This results in an expected redshift gradient
of $-10.9\times 10^{-19}/$cm, with the uncertainty bounded $<0.1\times
10^{-19}/$cm due to uncertainties in $g$ and $\Delta h$.

###  Normal operation, data taking, unblinding and analysis

We performed 14 blinded measurement runs of gravitational redshift data under
normal operating conditions over a 3-week campaign. Each run ranged in
duration from 1 to 4 hours, and was performed in conjunction with verification
of several experimental parameters to ensure that the associated systematic
effects are under control, such as the magnetic field gradient, density shift
coefficients, $\delta u$, and clock and lattice beam alignments (See Methods).
In each measurement run, the differential frequency of each ensemble pair was
extracted through ellipse fitting, with the associated Allan deviation
extrapolated to the full averaging time taken as the statistical uncertainty.
The corrections for density shifts and bias error from the ellipse fitting are
applied to each ensemble pair individually. The total uncertainty of each
clock comparison is given by the quadrature sum of its statistical uncertainty
and the uncertainties of systematic corrections. We analyze the measured
frequency gradient using two approaches.

In the first approach, the extracted frequency differences for 10 ensemble
pairs from each measurement run are plotted as a function of the height
differences. A linear fit is applied to each measurement run. The associated
fitted slopes from 14 measurement runs, after accounting for systematic
gradient corrections, are then weighted averaged, yielding a statistical
uncertainty of $0.7\times 10^{-19}/$cm, inflated by the square root of the
reduced $\chi^{2}$ statistic, $\chi^{2}_{\rm red}=1.62$ (Fig. 3a). Upon
completion of the measurements, the pseudo-randomly generated offset blinding
gradient was revealed and subtracted from the measurements. The blinding
offset gradient proved to be $+3.7\times 10^{-18}/$cm, and the measurements
before and after unblinding are shown in Extended Data Fig. 8. We find a
weighted mean frequency gradient of $[-12.4\pm 0.7_{(\rm stat)}\pm 2.4_{(\rm
sys)}]\times 10^{-19}/$cm, consistent with the expected redshift gradient of
$-10.9\times 10^{-19}/$cm within 1$\sigma$ total uncertainty.

In the second approach, we re-analyze the data and apply systematic
corrections to each pairwise clock comparison individually over the same raw
data set. For each pair, the total uncertainty is calculated as the quadrature
sum of the standard error of its weighted mean, systematic uncertainties that
don’t scale with height difference (density shift and ellipse fitting
corrections), and other systematic uncertainties that scale with height
difference. The weighted averaged frequency differences of each ensemble pair
are given in Extended Data Table 1. Through a final linear fit to the
differential frequencies as a function of the height differences, we find a
frequency gradient of $(-11.9\pm 2.5)\times 10^{-19}/$cm (Fig. 3b), again
fully consistent with the expected redshift gradient within $1\sigma$ total
uncertainty.

###  Clock height differences measurement with relativistic geodesy

Assuming that the gravitational redshift is given by Eq. 1, and treating the
ensemble array as a network of spatially distributed clocks with unknown
heights, we demonstrate relativistic geodesy below the geodetic limit [11, 12,
8] using synchronous clock comparisons (Extended Data Table 1). The height
difference $\Delta h$ for each clock pair is then given by

$\Delta h=\frac{\delta f}{f}\frac{c^{2}}{g},$ (21)

where $\delta f$ is the measured gravitational redshift, $g$ is the
independently measured local gravitational acceleration ($g=-9.803$ m$/$s2),
$f$ is the clock transition frequency, and $c$ is the speed of light. The
uncertainties of the extracted height differences are dominated by the
systematic uncertainties of the measured gravitational redshifts.

###  Data availability

The experimental data presented in this manuscript are available from the
corresponding author upon reasonable request.

###  Code availability

The code used for experimental control, data analysis, and simulation in this
work are available from the corresponding author upon reasonable request.

###  Acknowledgements

We acknowledge G.E.W. Marti, A. Jayich and T. Bothwell for fruitful
discussions and insightful comments on the manuscript. We acknowledge
technical contributions from B.N. Merriman, H. Li, V. Lochab and N. Ranabhat.
We are particularly grateful to B. Tikoff, E. M. Nelson and C. Ruggles from
the Department of Geoscience at the University of Wisconsin-Madison for
performing measurements of the gravitational acceleration in our laboratory.
This work was supported by the NIST Precision Measurement Grants program, the
Northwestern University Center for Fundamental Physics and the John Templeton
Foundation through a Fundamental Physics grant, the Wisconsin Alumni Research
Foundation, a Packard Fellowship for Science and Engineering, the Army
Research Office through agreement number W911NF-21-1-0012, and the National
Science Foundation under Grant No. 2143870.

###  Author contributions

All authors contributed to carrying out the experiments, data analysis and
writing the manuscript.

###  Competing interests

The authors declare no competing interests.

###  Correspondence and requests for materials

Correspondence and requests for materials should be addressed to Shimon
Kolkowitz.

## References

  * [1] Halir, R. & Flusser, J. Numerically stable direct least squares fitting of ellipses. Proc. of Sixth Intl Conf. Computer Graphics and Visualization 1 (1998)
  * [2] Estey, B. V. Ph. D. thesis (2016)
  * [3] Marti, G. E. et al. Imaging optical frequencies with 100 $\mu$Hz precision and 1.1 $\mu$m resolution. Phys. Rev. Lett. 120, 103201 (2008)
  * [4] Boyd, M. M. et al. Nuclear spin effects in optical lattice clocks. Phys. Rev. A 76, 022510 (2007)
  * [5] Boyd, M. M. Ph. D. thesis (2013)
  * [6] Beloy, K. et al. Atomic clock with $1\times 10^{-18}$ room-temperature blackbody stark uncertainty. Phys. Rev. Lett. 113, 260801 (2014)
  * [7] Ushijima, I., Takamoto, M., Das, M., Ohkubo, T., & Katori, H. Cryogenic optical lattice clocks. Nat. Photon. 9, 185–189 (2015)
  * [8] Beloy, K. et al. Faraday-shielded dc stark-shift-free optical lattice clock. Phys. Rev. Lett. 120, 183201 (2018)
  * [9] Lemonde, P. & Wolf, P. Optical lattice clock with atoms confined in a shallow trap. Phys. Rev. A 72, 033409 (2005)

Figure 1: $|$Atomic temperatures. a, Axial clock sideband thermometry taken at
130 $E_{\rm rec}$ lattice trap depth following in-lattice cooling. The reduced
height of the red sideband ($n_{z}\rightarrow n_{z}-1$) indicates that $99\%$
of atoms are populated in the lattice ground band ($n_{z}=0$). b, Radial
Doppler broadened profile taken at $15\leavevmode\nobreak\ {E_{\rm rec}}$. A
Gaussian fit (solid blue line) yields a radial temperature of 180 nK.

Figure 2: $|$Atomic coherence times. Decay of Ramsey contrast as a function of
dark times at $u_{op}=15\leavevmode\nobreak\ E_{\rm rec}$ (scatter points). An
exponential fit (solid line) yields an atomic coherence time of 32.1(9) s.

Figure 3: $|$Black body radiation shift evaluation. Evaluations of BBR shift
for all 10 ensemble pairs simultaneously. Solid lines are linear fits to the
scatter points. The y-axis scales are kept the same for all the plots. The
dense data points around 0.35 K correspond to temperature differences recorded
under normal operating conditions. The inset text list the ensemble pair
indices and corresponding height difference.

Figure 4: $|$Lattice light shift evaluation. a, Extraction of relative trap
depth difference $\delta u$ for 10 ensemble pairs. The extracted $\delta u$
are symmetric around 0, reflecting the Gaussian nature of the lattice beam
profile and that the ensembles are centered about the beam focus. b,
Evaluation of differential lattice light shifts at
$u_{op}=15\leavevmode\nobreak\ E_{\rm rec}$ for 10 ensemble pairs. The shifts
scale with spatial separation, but do not scale with $\delta u$, resulting in
a residual spatial lattice light shift gradient of $-8(1)\times
10^{-20}/E_{\rm rec}/$cm. c, Evaluation of the differential vector Stark
shift, which also doesn’t scale with $\delta u$. A linear fit (blue solid
line) yields a gradient of $-2.5(2)\times 10^{-18}/E_{\rm rec}/$cm.

Figure 5: $|$Probe Stark shift evaluation. a, Measured $\pi$ pulse durations
$T_{\pi}$ for ensemble 1 (red) and 5 (blue), the pair with largest spatial
separation (1 cm), as a function of the angle of clock beam misalignment. The
solid lines are quadratic fit to the data. b, Evaluated probe Stark shift at
the operational intensity $I_{op}$ at $\pm 500\leavevmode\nobreak\ \mu$rad
clock beam misalignment (estimated using camera). A linear fit to the data
places bounds on the probe Stark shift gradient at below $0.5\times
10^{-19}/$cm (the blue area represents a $\pm 50\leavevmode\nobreak\ \mu$rad
range).

Figure 6: $|$Ellipse fitting bias error evaluation. Monte-Carlo simulations of
bias error from ellipse fitting with input parameters of 2000 atoms per
ensemble, 85% contrast and 10 s Ramsey dark time (red solid line). The red
area represents our normal operating range of $(\pi/6,5\pi/6)$, where
fractional bias error is bounded to below $2\times 10^{-19}$ and is corrected
for each ellipse fitting during data analysis with uncertainty below
$0.5\times 10^{-19}$.

Figure 7: $|$Wannier-Stark ladder sidebands. a, Illustration of the Wannier-
Stark ladder in a 1D vertical optical lattice along the $\hat{z}$ direction.
The Wannier state at lattice site $n$ has energy $nmg\lambda_{L}/2$. The black
arrow represents the carrier clock transition where $n$ remains unchanged
($\Delta n=0$). The blue (red) arrow represents the first-order sideband
transition with $\Delta n=+1(-1)$. b, Representative plot of measured Wannier-
Stark sidebands at a shallow lattice depth ($5\leavevmode\nobreak\ E_{\rm
rec}$) recorded with Fourier-limited Rabi spectroscopy with a $\pi$ pulse
duration of 40 ms. The black line corresponds to the carrier ($\Delta n=0$)
transition, and the blue (red) line corresponds to off-site transitions,
$\Delta n=+1(-1)$. The splitting between $\Delta n=\pm 1$ transitions is used
to extract the lattice tilt.

Figure 8: $|$Measurements before and after removal of the blinded offset
gradient. a, Measured fractional frequency gradients (blue scatter points)
with blinded offset. The blue area represents $\pm 1\sigma$ total uncertainty.
The black solid line is the expected gravitational redshift gradient of
$-10.9\times 10^{-19}/$cm. The green area represents the range of possible
outcomes upon removal of the blinded offset gradient, corresponding to a $\pm
5\times 10^{-18}/$cm range. All corrections and uncertainties were finalized
at this stage, without knowledge of the blinded offset. b, Measured fractional
frequency gradients following removal of the blinded offset. Red scatter
points are the same data set as in a with the blinded offset gradient
($+3.7\times 10^{-18}/$cm) removed. This represents the same data shown in
Fig. 3a, but here the scale of the y-axis is kept the same as in panel a for
the sake of comparison.

Figure 1: $|$Clock height difference measurement. The measured gravitational
redshift of each clock comparison is used to extract the height difference
given the independently measured local gravitational acceleration, $g=-9.803$
m$/$s2.