# All-in-one Quantum Diamond Microscope for Sensor Characterization.

Connor Roncaioli<EMAIL_ADDRESS>DEVCOM Army Research
Laboratory  Quantum Technology Center, University of Maryland  Connor Hart
<EMAIL_ADDRESS>Quantum Technology Center, University of Maryland  Ronald
Walsworth<EMAIL_ADDRESS>Quantum Technology Center, University of Maryland
Donald P. Fahey<EMAIL_ADDRESS>DEVCOM Army Research Laboratory
Quantum Technology Center, University of Maryland

###### Abstract

Nitrogen-vacancy (NV) centers in diamond are a leading modality for magnetic
sensing and imaging under ambient conditions. However, these sensors suffer
from degraded performance due to paramagnetic impurities and regions of stress
in the diamond crystal lattice. This work demonstrates a quantum diamond
microscope (QDM) for the simultaneous mapping and spatial correlation of key
properties of a millimeter-scale NV-diamond sensor chip, including: NV
ensemble photoluminescence (PL) amplitude, spin-lattice relaxation time (T1),
and homogeneous and inhomogeneous spin coherence lifetimes (T2 and
T${}_{2}^{*}$), as well as lattice stress/strain, birefringence magnitude, and
birefringence angle of the diamond crystal.

††preprint: AIP/123-QED

##  I Introduction

Solid state quantum sensors, such as the nitrogen-vacancy (NV) center in
diamond, are undergoing rapid research and commercial expansion. Example NV-
diamond measurement modalities include compact scalarTaylor  _et al._ (2008);
Rondin  _et al._ (2014); Eisenach  _et al._ (2021); Fescenko  _et al._ (2020);
Ziem  _et al._ (2019); Aslam  _et al._ (2023) and vectorSchloss  _et al._
(2018); Wang  _et al._ (2015) magnetic field sensing and imaging, stable clock
oscillatorsHodges  _et al._ (2013); Trushein  _et al._ , electric field
sensorsIwasaki  _et al._ (2017); Scholten  _et al._ (2021), nanodiamond
temperature sensorsNeumann  _et al._ (2013); Aslam  _et al._ (2023), mapping
of strain in the diamond crystalKehayias  _et al._ (2019); Barfuss  _et al._
(2019), microwave mode coolingZhang  _et al._ (2022); Fahey  _et al._ ,
extreme environment magnetometryToyli  _et al._ (2012), and pressure
sensing.Bhattacharyya  _et al._ (2024)

A widely used platform is the quantum diamond microscope (QDM), which employs
a dense surface layer of NVs on a macroscopic diamond plate to enable wide-
field micron-scale vector imaging of magnetic fields at very low standoff
distances.Levine  _et al._ (2019); Barry  _et al._ (2020) Useful operation of
the QDM requires detailed characterization of NV-diamond sensor properties for
particular use casesAlghannam and Hemmer (2019); Osterkamp  _et al._ (2019);
Ofori-Okai  _et al._ (2012), as well as to identify imperfections that can
limit QDM performance.Kehayias  _et al._ (2019); Marshall  _et al._ (2022);
Knauer, Hadden, and Rarity (2020) Here we demonstrate an "all-in-one QDM" that
can image several parameters relevant to NV ensemble sensing across a
millimeter diamond plate, including NV spin resonance lineshift, decoherence,
depolarization time scales, stress/strain, and birefringence. Measurements are
co-registered to simplify analysis of complex features that appear across
modalities.

Figure 1:  Schematic of all-in-one QDM for NV-diamond sensor characterization.
Two optical excitation pathways are exchanged by use of a flipper mirror, with
an epifluorescent option and a non-focused full-sample excitation option.
Microwaves (MWs) are coupled in to an edge inductive gap of the loop gap
resonator (LGR), creating a strong, homogeneous MW field across the full NV-
diamond sample. Permanent magnets create a uniform bias magnetic field that
Zeeman shifts NV spin resonances in order to isolate or overlap resonance
frequencies as needed. Optical objectives above and below the sample are used
for epifluorescent or full field-of-view measurements. The avalanche
photodiode (APD) and optical spectrometer are coupled to the above-sample
collection pathway. The APD is used for time-resolved NV photoluminescence
(PL) measurements, while the spectrometer is used for MW assisted NV charge
state measurements. The below-sample objective is coupled to a long-pass
filter and camera used to image PL for spatially resolved NV characterization.
An LED with linear polarization control and appropriate circular polarization
filter can be switched for the top objective and bottom filters, respectively,
to create a transmission measurement of optical birefringence using the same
below-sample objective and camera.

##  II Instrument Overview

The all-in-one QDM is configurable between multiple measurement modalities, as
listed in Table 1. A high-level depiction of the apparatus is illustrated in
Fig. 1 and the supplemental material includes a summary of each measurement
modality and list of components for each subsystem.

Laser light for NV excitation can be sent through a top objective for point-
excitation, or to a periscope that illuminates a mm-scale region of the
diamond along a shallow angle. Similarly there are two collection pathways for
NV photoluminescence (PL), the top objective can be used to collect
epifluorescence and spatially filter it before being sent to either a
photodiode for time resolved confocal measurement or a fiber-coupled
spectrometer for spectroscopy; or the PL can be collected by an imaging
objective mounted beneath the sample space, which directs light to a camera
for wide-field imaging. This combination of options affords great flexibility:
e.g., the all-in-one QDM can collect and correlate confocal and imaged NV
measurements. The top objective can be exchanged with a
$530\text{\,}\mathrm{nm}$ LED with 360°linear polarization control, and the
filters in the imaging path can be switched with a circular polarization
analyzer to change the system into an optical birefringence transmission
measurement, as detailed in Section III.3.

The system is capable of strong MW driving of NV spin transitions
(>$1\text{\,}\mathrm{MHz}$ Rabi frequency) across the full field of view of
the imaging objective ($\approx$$5\text{\,}\mathrm{mm}$). Achieving an NV Rabi
frequency $\approx\frac{1}{\textrm{T}_{2}^{*}}$ across the entire diamond
sample is necessary to enable simultaneous imaging of key properties such as
the NV ensemble spectrum of optically detected magnetic resonance (ODMR). This
MW capability is enabled by the use of a metallic loop gap resonator (LGR)
innovated at MIT-LL.Eisenach  _et al._ (2018) Details are provided in the
supplemental material.

Figure 2:  a) Example map of NV ensemble Rabi frequency produced by all-in-one QDM across $\approx$3x$3\text{\,}\mathrm{mm}$ region of a diamond sensor plate. NV-diamond sample is described in section IV. Observed variation in Rabi frequency across field-of-view results from spatial gradients in MW coupling strength to NV spin transitions, e.g., due to inhomogeneous MW drive amplitude or angle of incidence relative to the interrogated NV dipole, or any cause of NV spin resonance lineshift or broadening. Spatial scale bar is shown. b) Histogram of fit values provides an estimate of Rabi frequency distribution across all-in-one QDM field-of-view. 10th and 90th percentiles, marked with dashed black lines, fall within $\pm$10% of mean Rabi frequency value. Inset: example single pixel Rabi oscillation data set and associated fit in blue.  Table 1: All-in-one QDM measurement modalities: CLE - confocal laser excitation using epifluorescent illumination, FSLE - full-sample laser excitation, LED - light-emitting diode, PD - photodiode, CMOS - complementary metal-oxide semiconductor camera. | Measurement Details |
---|---|---
Measurement | Pump Geometry | Read-Out Configuration | Purpose
Constant wave-optically detected | CLE | PD | MW resonance frequencies
magnetic resonance (CW-ODMR) |  |  | of NV spin transitions
Rabi | CLE | PD | MW driving rate of NV spins
T1 Relaxometry | CLE | PD | NV spin-lattice relaxation time (T1)
Ramsey | CLE | PD | NV spin dephasing time (T${}_{2}^{*}$)
Hahn echo | CLE | PD | NV spin decoherence time (T2)
Dynamic decoupling | CLE | PD | Performance of noise-cancelling
|  |  | pulse sequences outside of
|  |  | specified frequency window
Photoluminescence (PL) distribution | FSLE or CLE | CMOS | NV PL intensity and epi-spot location111A supplement to photodiode measurements, showing the precise location of the pump spot relative to the sample in imaged measurements
Imaged CW-ODMR | FSLE | CMOS | NV spin resonance across sample
Imaged Rabi | FSLE | CMOS | MW driving rate of NVs across sample
Imaged T1 relaxometry | FSLE | CMOS | NV T1 across sample
Imaged Ramsey | FSLE | CMOS | NV T${}_{2}^{*}$ across sample
Imaged Hahn echo | FSLE | CMOS | NV T2 across sample
Imaged dynamic decoupling | FSLE | CMOS | Performance of noise-cancelling
|  |  | pulse sequences across sample
| Illumination Source |  |
Imaged quantitative birefringence | LED | CMOS | Spatial distribution of
|  |  | stress/strain within diamond plate

In order to produce spatially resolved maps of parameters of interest, with
micron-scale resolution across a mm-scale field-of-view, upwards of 1 million
distinct data sets require fitting. This analysis requires the use of parallel
GPU-based fitting software, such as GPUfit and the associated python-based
pyGPUfit module.Przybylski  _et al._ (2017) For CW-ODMR, Rabi, and T1
relaxometry, measurements are reported as contrast:

$\textrm{contrast}=\frac{\textrm{I}_{\textrm{MWon}}}{\textrm{I}_{\textrm{MWoff}}}$
(1)

Here, I${}_{\textrm{MWon}}$ is the measured PL signal intensity with MWs on
and I${}_{\textrm{MWoff}}$ is the reference PL intensity measurement with MWs
off. For phase-sensitive modalities (Ramsey and Hahn echo) values are reported
as visibility:

$\textrm{visibility}=\frac{\textrm{I}_{\textrm{MW+}}-\textrm{I}_{\textrm{MW-}}}{\textrm{I}_{\textrm{MW+}}+\textrm{I}_{\textrm{MW-}}}$
(2)

Here, I${}_{\textrm{MW}+}$ is the measured PL intensity where all MW pulses
applied to the NVs have the same phase, and I${}_{\textrm{MW}-}$ is PL
intensity with the final $\frac{\pi}{2}$ MW pulse of opposite phase. More
details, such as the functional form of the various fit procedures and
techniques for evaluating fit appropriateness are provided in the supplemental
material.

##  III Experimental Modalities

Figure 3: Example NV ensemble measurements for diamond sample described in
section IV, and associated fits for data, collected with all-in-one QDM in
epifluorescent configuration. a) CW-ODMR spectra; b) Rabi oscillation; c)
Ramsey measurement; d) Hahn echo; e) T1 relaxation; and f) NV PL spectra
demonstrating MW-assisted charge state decomposition.Aude Craik  _et al._
(2020) For the CW-ODMR data, LGR resonant frequency $\approx$
$2785\text{\,}\mathrm{MHz}$ and a bandwidth $\approx$
$80\text{\,}\mathrm{MHz}$, which enhances MW driving of nearby resonant NV
spectral features arising from different projections of bias magnetic field on
four NV orientations within the diamond crystal.

###  III.1 CW-ODMR

Confocal continuous wave (CW)-ODMR is used for high-sensitivity measurement of
NV spin resonance frequencies at a near-diffraction-limited point and focal
depth. Imaged CW-ODMR, performed via non-focused full-sample laser
illumination, allows measurement of the NV spin resonance spectrum on a pixel-
by-pixel basis, but necessarily collects NV PL from all depths of field. This
modality can map interactions that affect the NV spin Hamiltonian, including
stress within the diamond. Based on previous work Kehayias  _et al._ (2019);
Barfuss  _et al._ (2019) we describe the NV stress tensor:

$\displaystyle\sigma_{diag}=\frac{1}{4a_{1}}[M_{z,1}+M_{z,2}+M_{z,3}+M_{z,4}]$
(3)
$\displaystyle\sigma_{XY}=\frac{1}{8a_{2}}[M_{z,1}+M_{z,2}-M_{z,3}-M_{z,4}]$
(4)
$\displaystyle\sigma_{XZ}=\frac{1}{8a_{2}}[M_{z,1}-M_{z,2}+M_{z,3}-M_{z,4}]$
(5)
$\displaystyle\sigma_{YZ}=\frac{1}{8a_{2}}[M_{z,1}-M_{z,2}-M_{z,3}+M_{z,4}]$
(6)

where $\sigma_{diag}=\sigma_{XX}+\sigma_{YY}+\sigma_{ZZ}$ is the diagonal
stress tensor component, the individual values of which are obscured in
present QDM measurements; Mz,i are the asymmetric lineshift values for each NV
orientation classes within the diamond crystal about their mean value; and
$\\{\textrm{a}_{1},\textrm{a}_{2}\\}=\\{\textrm{4.86},\textrm{-3.7}\\}$
$\frac{$\text{\,}\mathrm{MHz}$}{$\text{\,}\mathrm{GPa}$}$ is the spin-stress
coupling constant. Using these relationships we can partially reconstruct the
stress tensor, including the shear components, from QDM measurements (see
supplemental material). Note that determination of the individual uniaxial
terms ($\sigma_{XX}$, etc.) is possible by changing the bias magnetic field
direction, though this is not implemented here.Kehayias  _et al._ (2019)

Figure 4:  Example of NV-diamond sample characterization performed with all-
in-one QDM and displayed as image data. a) $\sigma_{diag}$ extracted from CW-
ODMR; b) T${}_{2}^{*}$ from Ramsey; c) T1 from relaxometry; d) T2 from Hahn
echo; e) stress magnitude and f) stress angle determined from optical
birefringence. Correlations in a)-f) relate to regions of high stress and
stress gradients.  Figure 5:  Expanded-view QDM characterization images for
region of NV-diamond sample indicated by box in Fig. 4a): this region of
enhanced stress appears as a channel between two highly localized linear
regions. a) Diagonal stress tensor component $\sigma_{diag}$ extracted from
CW-ODMR corresponds to isometric stress. b-d) Off diagonal stress tensor
components correspond to shear stress, with $\sigma_{XY}$ being in-plane shear
and $\sigma_{YZ}$ and $\sigma_{ZX}$ being out-of-plane shear components. e)
$\textrm{T}_{2}^{*}$ map shows qualitative correlation with stress
measurements. f) $\textrm{T}_{2}$ map is near-uniform with the exception of a
few marked regions. g) stress magnitude and h) stress angle determined from
optical birefringence measurements.

###  III.2 Pulsed Protocols

Pulsed NV measurements are perfomed in either the confocal or wide-field
imaging modes, including protocols for Rabi oscillation, Ramsey, Hahn echo,
and $\textrm{T}_{1}$ relaxometry, see Figs. 3 and 4. Importantly, the all-in-
one QDM allows such measurements over a millimeter-scale field-of-view —
whereas past measurements were limited to
$\approx$$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$Levine  _et al._ (2019) —
due to integration of the LGR to provide strong, homogeneous MW driving of NV
spin transitions.

Imaging the NV Rabi frequency allows evaluation of MW field homogeneity.
Ramsey and Hahn echo measurements are used to evaluate the NV ensemble spin
dephasing time (T${}_{2}^{*}$); and are conducted at a fixed MW detuning
($\approx$ $5\text{\,}\mathrm{MHz}$) from the NV center frequency in order to
have sufficiently fast oscillation to properly fit the decay envelope. Hahn
echo measurements are used to measure the NV ensemble decoherence time (T2).
Both Ramsey and Hahn echo measurements are collected with phase dependence of
the final MW pulse in order to measure visibility instead of contrast for
improved noise cancellation. Finally T1 relaxometry provides an upper limit to
NV spin lifetimes in the sample. Dynamical decoupling schemes are feasible for
this system but are not implemented in the present work. Examples of imaged
fitting parameters extracted from these measurements are provided in Figs. 4
and 5. Correlation is observed between high-contrast features in the maps of
stress and T${}_{2}^{*}$, which is consistent with past studies of CW-ODMR and
Ramsey microscopy of NV-diamond samples.Marshall  _et al._ (2022) However, the
T2 and T1 maps are largely homogeneous, as anticipated for the high-quality
NV-diamond sample studied here. More detail is provided in section IV.

###  III.3  Optical Birefringence and Stress Imaging

The cubic symmetry of ideal diamond has zero optical birefringence; however,
real diamond samples typically have local stress that breaks this symmetry,
resulting in observable birefringence.Lang (1967) The all-in-one QDM provides
a quantitative measure of stress in the diamond sample via a switchable stage
that brings an LED (wavelength = $\sigma$ = $530\text{\,}\mathrm{nm}$) with a
rotating linear polarizer above the diamond sample and a circular analyzer
after the below-sample objective. In this operational mode, optical
transmission measurements yield the full 2D anisotropic retardance vector
without adjusting the sample.Kehayias  _et al._ (2019) The stress angle $\phi$
and retardance phase $\delta=\frac{2\pi}{\lambda}\Delta\textrm{nL}$, where L
is the sample thickness and $\Delta n$ is the difference in refractive index
between the two polarizations, are related to the optical transmission
intensity by the following equation:

$\textrm{I}_{\textrm{i}}=\frac{1}{2}\textrm{I}_{\textrm{0}}[1+\textrm{sin}2(\alpha_{\textrm{i}}-\phi)\textrm{sin}\delta]$
(7)

By measuring I${}_{\textrm{i}}$ as the linear polarizer angle
$\alpha_{\textrm{i}}$ rotates through 180°, $\phi$ and sin $\delta$ are
determined. $\delta$ is related to the optical wavelength $\lambda$ and the
diamond stress magnitude (|$\sigma$|) and the index of refraction n=2.42 by
the equation:

$|\sigma|\approx\frac{2}{3\pi}\frac{|\delta|\lambda}{\textrm{Ln}^{3}\textrm{q}_{\textrm{iso}}}$
(8)

where q${}_{\textrm{iso}}$=0.3x10-12$\text{\,}\mathrm{Pa}$-1 is a piezo-
optical coefficient.Howell  _et al._ (2010)

Stress determined in this way may be buried deep in the undoped diamond, far
away from the active layer of NVs. Additionally, small regions of intense
stress in the NV layer may not be detectable by birefringence against the
contribution of a large undoped diamond, if the active layer is comparatively
small. Furthermore, if the diamond has significant strain or the sample is of
sufficient thickness, the optical birefringence measurement can become phase-
ambiguous once the difference between the slow and fast axes
$|\delta|\geq\frac{\lambda}{2}$. These are all reasons that multi-modal
detection is crucial for complete characterization of NV samples. Comparison
of measurements of stress determined by CW-ODMR (contribution confined to the
NV layer) to optical birefringence (contribution from the entire diamond) are
presented in Figs. 4a), 4e), and 4f).

###  III.4 Optical Spectroscopy and NV Charge State Characterization

The all-in-one QDM can perform optical spectroscopy from a diamond sample
excited using either confocal or shallow-angle excitation. Modulating NV PL
spectra using resonant MWs allows determination of the relative concentration
of NV charge states (NV${}^{\textrm{-}}$ and NV${}^{\textrm{0}}$).Aude Craik
_et al._ (2020) Note that only the NV${}^{\textrm{-}}$ charge state exhibits
spin-state-dependent PL, and hence is useful for quantum sensing.

NV charge states are metastable and dependent on both intrinsic material and
extrinsic experimental conditions. Intrinsic conditions that can affect NV
charge state include local availability of electron donors by dopingHerbschleb
_et al._ (2019) and local chemical potentialKaraveli  _et al._ (2016). Some
samples can have coexisting regions of different charge state ratios.Savvin
_et al._ (2021) Extrinsic conditions that can affect the NV charge state
include electric field bias,Lozovoi  _et al._ (2020)Lozovoi  _et al._ (2023)
laser intensity,Manson and Harrison (2005) and excitation wavelength.Aslam
_et al._ (2013) NV${}^{\textrm{-}}$ fraction as a function of laser intensity
is particularly important to characterize since it falls off at high
intensities, which needs to be balanced against NV PL for optimal
sensing.Barry  _et al._ (2020)

##  IV Diamond Case Study

As a demonstration of the all-in-one QDM, we characterize a quantum grade NV-
diamond sample produced by Element Six. This sample is 3x3x0.5mm, with a
$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$ thick 16ppm N, 3ppm NV layer. The
sample is treated with a tri-acid perchlorate solutionBrown  _et al._ (2019)
in order to remove surface contaminants and characterized by XPS to discern
surface composition; additional details are provided in the supplemental
material. The sample is then characterized with the all-in-one QDM, performing
both confocal and imaged CW-ODMR, Rabi, Ramsey, and Hahn echo, as well as
imaged optical birefringence of the entire sample, imaged stress of the NV
layer, and confocal MW-assisted spectroscopy to determine relative NV charge
states (see Figs. 2, 3, 4 and 5). After alignment, the diamond sample remains
untouched until all measurements are completed, so that imaged modalities are
coregistered to within one pixel, enabling precise correlation of features
between measurements.

For the wide-field imaging modalities, significant heating of the diamond
sample is caused by laser illumination at typical input power of a few watts.
Heating leads to an asymmetric lineshift of NV spin resonances, which can
result in a lineshift gradient across the sample. To mitigate the effect of
the heat load, the diamond sample is mounted cantilever style to a thermally
conductive SiC plate; and a $5\text{\,}\mathrm{\SIUnitSymbolMicro s}$ delay is
placed between the laser initialization pulse and the MW protocol used for NV
measurements. Cantilever mounting also simplifies the subtraction of the
birefringence background, since the only contribution comes from free space
and rigidly aligned optics.

Fig. 4 presents example images of the diamond samples’s NV ensemble asymmetric
resonance lineshift in CW-ODMR (corresponding to the isometric stress
component $\sigma_{diag}$ within the NV layer), T${}_{2}^{*}$, T1, and T2; as
well as optical birefringence over the full sample thickness, yielding both
the relative stress magnitude and angle. Regions of stress typically
correspond to NV spin resonance broadening and reduced T2, and T${}_{2}^{*}$,
limiting the sensitivity to magnetic fields in these locations.Kehayias  _et
al._ (2019); Barfuss  _et al._ (2019) Over 90% of the sample has measured
T${}_{2}^{*}$ within 10% of the median value; and almost 99% of the sample has
T2 within 5% of the median T2 value. Avoiding the few regions of enhanced
stress and reduced T${}_{2}^{*}$ and T2 is thereby enabled by all-in-on QDM
characterization.

Figure 5 shows a specific region of enhanced stress at higher spatial
resolution, selected from the QDM images of Fig. 4. In Fig. 5a) the isometric
stress component ($\sigma_{diag}$) shows a region of increased stress that is
between two linear regions in a channel
$\approx$$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ across. Additional stress
petals are visible nearby as much smaller features that are tens of microns
across. Fig. 5b) shows that the two edges of the region of enhanced stress
have related, but opposite, polarity shear components that extend across the
region of increased $\sigma_{diag}$. Figs. 5c) and 5d) demonstrate that the
out-of-plane shear components are localized only to the edges of the region of
increased stress, and do not extend into the internal region. Figs. 5e) and
5f) show that increased isometric stress is correlated with (and likely causal
to) significantly reduced T${}_{2}^{*}$, but with minimal effect on T2. In
particular, the edges of the stress region, which correspond to localized
behavior in the shear components or to steep gradients in isometric stress,
correlate with decreased T${}_{2}^{*}$ values. Figs. 5g) and 5h) show stress
magnitude and phase, respectively, determined from optical birefringence
measurements; and provide quanitative confirmation of high stress regions
determined from NV measurements.

##  V Conclusion and Future Work

We present a multi-functional quantum diamond microscope (QDM) for
comprehensive characterization of NV-diamond samples. The apparatus employs
automated components that enable remote operation. With modest future
development, we envision push-button operation of all modalities described in
this paper without user intervention. More significant future effort will
likely be needed to extend the capability of this "all-in-one QDM" to remote
operation on multiple samples in sequence, as a result of challenges with
careful sample alignment, background subtraction with regards to
birefringence, as well as good thermal heat-sinking.

The apparatus described in this paper should be extendable to any optically
driven and detected solid-state quantum defect that operates at room
temperature. There is a vast library of potential defects to be explored for
optically driven quantum sensing, including vacancies and site defectsFalk
_et al._ (2013), transition metalDiler  _et al._ (2020) and rare earth
dopantsZhang  _et al._ (2020)Spindlberger  _et al._ (2019) in both diamond and
other semiconductor materials. The scale of characterization demands will
continue to grow alongside developments in semiconductor defect systems,
underlining the need for rapid characterization systems as described in this
paper.

###### Acknowledgements.

The authors thank Dr. Brenda VanMil, Dr. Sean Blakley and Dr. Matthew Trusheim
for insightful conversations on solid state systems and NV physics. We also
thank Dr. Adam Biacchi at NIST for implementing the conventional fumehood safe
tri-acid diamond cleaning method used in this paper, as well as XPS analysis
of the sample surface to confirm a successful cleaning and to interrogate
surface chemistry. Support for this work was provided by the U.S. Army
Research Laboratory, including under Contract No. W911NF1920181, and the
University of Maryland Quantum Technology Center.

## References

  * Taylor  _et al._ (2008) J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin, “High-sensitivity diamond magnetometer with nanoscale resolution,” Nature Physics 4, 810–816 (2008).
  * Rondin  _et al._ (2014) L. Rondin, J.-P. Tetienne, T. Hingant, J.-F. Roch, P. Maletinsky, and V. Jacques, “Magnetometry with nitrogen-vacancy defects in diamond,” Reports on Progress in Physics 77, 056503 (2014).
  * Eisenach  _et al._ (2021) E. R. Eisenach, J. F. Barry, M. F. O’Keeffe, J. M. Schloss, M. H. Steinecker, D. R. Englund, and D. A. Braje, “Cavity-enhanced microwave readout of a solid-state spin sensor,” Nature Communications 12 (2021), 10.1038/s41467-021-21256-7.
  * Fescenko  _et al._ (2020) I. Fescenko, A. Jarmola, I. Savukov, P. Kehayias, J. Smits, J. Damron, N. Ristoff, N. Mosavian, and V. M. Acosta, “Diamond magnetometer enhanced by ferrite flux concentrators,” Physical Review Research 2 (2020), 10.1103/PhysRevResearch.2.023394.
  * Ziem  _et al._ (2019) F. Ziem, M. Garsi, H. Fedder, and J. Wrachtrup, “Quantitative nanoscale mri with a wide field of view,” Scientific Reports 9 (2019), 10.1038/s41598-019-47084-w.
  * Aslam  _et al._ (2023) N. Aslam, H. Zhou, E. K. Urbach, M. J. Turner, R. L. Walsworth, M. D. Lukin, and H. Park, “Quantum sensors for biomedical applications,” Nature Reviews Physics 5, 157–169 (2023).
  * Schloss  _et al._ (2018) J. M. Schloss, J. F. Barry, M. J. Turner, and R. L. Walsworth, “Simultaneous broadband vector magnetometry using solid-state spins,” Physical Review Applied 10 (2018), 10.1103/PhysRevApplied.10.034044.
  * Wang  _et al._ (2015) P. Wang, Z. Yuan, P. Huang, X. Rong, M. Wang, X. Xu, C. Duan, C. Ju, F. Shi, and J. Du, “High-resolution vector microwave magnetometry based on solid-state spins in diamond,” Nature Communications 6, 6631 (2015).
  * Hodges  _et al._ (2013) J. S. Hodges, N. Y. Yao, D. Maclaurin, C. Rastogi, M. D. Lukin, and D. Englund, “Timekeeping with electron spin states in diamond,” Phys. Rev. A 87, 032118 (2013).
  * (10) M. E. Trushein, K. Jacobs, J. E. Hoffman, D. P. Fahey, D. A. Braje, and D. Englund, “A polariton-stabilized spin clock,” arXiv 10.48550/arXiv.2009.02427.
  * Iwasaki  _et al._ (2017) T. Iwasaki, W. Naruki, K. Tahara, T. Makino, H. Kato, M. Ogura, D. Takeuchi, S. Yamasaki, and M. Hatano, “Direct nanoscale sensing of the internal electric field in operating semiconductor devices using single electron spins,” ACS Nano 11, 1238–1245 (2017).
  * Scholten  _et al._ (2021) S. C. Scholten, A. J. Healey, I. O. Robertson, G. J. Abrahams, D. A. Broadway, and J.-P. Tetienne, “Widefield quantum microscopy with nitrogen-vacancy centers in diamond: Strengths, limitations, and prospects,” Journal of Applied Physics 130, 150902 (2021), https://doi.org/10.1063/5.0066733 .
  * Neumann  _et al._ (2013) P. Neumann, I. Jakobi, F. Dolde, C. Burk, R. Reuter, G. Waldherr, J. Honert, T. Wolf, A. Brunner, J. H. Shim, D. Suter, H. Sumiya, J. Isoya, and J. Wrachtrup, “High-Precision nanoscale temperature sensing using single defects in diamond,” Nano Lett. 13, 2738–2742 (2013).
  * Kehayias  _et al._ (2019) P. Kehayias, M. J. Turner, R. Trubko, J. M. Schloss, C. A. Hart, M. Wesson, D. R. Glenn, and R. L. Walsworth, “Imaging crystal stress in diamond using ensembles of nitrogen-vacancy centers,” Phys. Rev. B 100, 174103 (2019).
  * Barfuss  _et al._ (2019) A. Barfuss, M. Kasperczyk, J. Kölbl, and P. Maletinsky, “Spin-stress and spin-strain coupling in diamond-based hybrid spin oscillator systems,” Phys. Rev. B 99, 174102 (2019).
  * Zhang  _et al._ (2022) Y. Zhang, Q. Wu, H. Wu, X. Yang, S.-L. Su, C. Shan, and K. Mølmer, “Microwave mode cooling and cavity quantum electrodynamics effects at room temperature with optically cooled nitrogen-vacancy center spins,” npj Quantum Information 8, 125 (2022).
  * (17) D. P. Fahey, K. Jacobs, M. J. Turner, H. Choi, J. E. Hoffman, D. Englund, and M. E. Trusheim, “Steady-state microwave mode cooling with a diamond nv ensemble,” arXiv 10.48550/arXiv.2203.03462.
  * Toyli  _et al._ (2012) D. M. Toyli, D. J. Christle, A. Alkauskas, B. B. Buckley, C. G. V. de Walle, and D. D. Awschalom, “Measurement and control of single nitrogen-vacancy center spins above 600 k,” Physical Review X 2 (2012), 10.1103/PhysRevX.2.031001.
  * Bhattacharyya  _et al._ (2024) P. Bhattacharyya, W. Chen, X. Huang, S. Chatterjee, B. Huang, B. Kobrin, Y. Lyu, T. J. Smart, M. Block, E. Wang, Z. Wang, W. Wu, S. Hsieh, H. Ma, S. Mandyam, B. Chen, E. Davis, Z. M. Geballe, C. Zu, V. Struzhkin, R. Jeanloz, J. E. Moore, T. Cui, G. Galli, B. I. Halperin, C. R. Laumann, and N. Y. Yao, “Imaging the meissner effect in hydride superconductors using quantum sensors,” Nature 627, 73–79 (2024).
  * Levine  _et al._ (2019) E. V. Levine, M. J. Turner, P. Kehayias, C. A. Hart, N. Langellier, R. Trubko, D. R. Glenn, R. R. Fu, and R. L. Walsworth, “Principles and techniques of the quantum diamond microscope,” Nanophotonics 8, 1945–1973 (2019).
  * Barry  _et al._ (2020) J. F. Barry, J. M. Schloss, E. Bauch, M. J. Turner, C. A. Hart, L. M. Pham, and R. L. Walsworth, “Sensitivity optimization for nv-diamond magnetometry,” Rev. Mod. Phys. 92, 015004 (2020).
  * Alghannam and Hemmer (2019) F. Alghannam and P. Hemmer, “Engineering of shallow layers of nitrogen vacancy colour centres in diamond using plasma immersion ion implantation,” Scientific Reports 9, 5870 (2019).
  * Osterkamp  _et al._ (2019) C. Osterkamp, M. Mangold, J. Lang, P. Balasubramanian, T. Teraji, B. Naydenov, and F. Jelezko, “Engineering preferentially-aligned nitrogen-vacancy centre ensembles in cvd grown diamond,” Scientific Reports 9, 5786 (2019).
  * Ofori-Okai  _et al._ (2012) B. K. Ofori-Okai, S. Pezzagna, K. Chang, M. Loretz, R. Schirhagl, Y. Tao, B. A. Moores, K. Groot-Berning, J. Meijer, and C. L. Degen, “Spin properties of very shallow nitrogen vacancy defects in diamond,” Phys. Rev. B 86, 081406 (2012).
  * Marshall  _et al._ (2022) M. C. Marshall, R. Ebadi, C. Hart, M. J. Turner, M. J. Ku, D. F. Phillips, and R. L. Walsworth, “High-precision mapping of diamond crystal strain using quantum interferometry,” Phys. Rev. Appl. 17, 024041 (2022).
  * Knauer, Hadden, and Rarity (2020) S. Knauer, J. P. Hadden, and J. G. Rarity, “In-situ measurements of fabrication induced strain in diamond photonic-structures using intrinsic colour centres,” npj Quantum Information 6, 50 (2020).
  * Eisenach  _et al._ (2018) E. R. Eisenach, J. F. Barry, L. M. Pham, R. G. Rojas, D. R. Englund, and D. A. Braje, “Broadband loop gap resonator for nitrogen vacancy centers in diamond,” Review of Scientific Instruments 89 (2018), 10.1063/1.5037465.
  * Przybylski  _et al._ (2017) A. Przybylski, B. Thiel, J. Keller-Findeisen, B. Stock, and M. Bates, “Gpufit: An open-source toolkit for gpu-accelerated curve fitting,” Scientific Reports 7, 15722 (2017).
  * Aude Craik  _et al._ (2020) D. Aude Craik, P. Kehayias, A. Greenspon, X. Zhang, M. Turner, J. Schloss, E. Bauch, C. Hart, E. Hu, and R. Walsworth, “Microwave-assisted spectroscopy technique for studying charge state in nitrogen-vacancy ensembles in diamond,” Phys. Rev. Applied 14, 014009 (2020).
  * Lang (1967) A. R. Lang, “Causes of birefringence in diamond,” Nature 213, 248–251 (1967).
  * Howell  _et al._ (2010) D. Howell, I. G. Wood, D. P. Dobson, A. P. Jones, L. Nasdala, and J. W. Harris, “Quantifying strain birefringence halos around inclusions in diamond,” Contributions to Mineralogy and Petrology 160, 705–717 (2010).
  * Herbschleb  _et al._ (2019) E. D. Herbschleb, H. Kato, Y. Maruyama, T. Danjo, T. Makino, S. Yamasaki, I. Ohki, K. Hayashi, H. Morishita, M. Fujiwara, and N. Mizuochi, “Ultra-long coherence times amongst room-temperature solid-state spins,” Nature Communications 10, 3766 (2019).
  * Karaveli  _et al._ (2016) S. Karaveli, O. Gaathon, A. Wolcott, R. Sakakibara, O. A. Shemesh, D. S. Peterka, E. S. Boyden, J. S. Owen, R. Yuste, and D. Englund, “Modulation of nitrogen vacancy charge state and fluorescence in nanodiamonds using electrochemical potential,” Proceedings of the National Academy of Sciences 113, 3938–3943 (2016), https://www.pnas.org/doi/pdf/10.1073/pnas.1504451113 .
  * Savvin  _et al._ (2021) A. Savvin, A. Dormidonov, E. Smetanina, V. Mitrokhin, E. Lipatov, D. Genin, S. Potanin, A. Yelisseyev, and V. Vins, “Nv– diamond laser,” Nature Communications 12, 7118 (2021).
  * Lozovoi  _et al._ (2020) A. Lozovoi, H. Jayakumar, D. Daw, A. Lakra, and C. A. Meriles, “Probing metastable space-charge potentials in a wide band gap semiconductor,” Phys. Rev. Lett. 125, 256602 (2020).
  * Lozovoi  _et al._ (2023) A. Lozovoi, Y. Chen, G. Vizkelethy, E. Bielejec, J. Flick, M. W. Doherty, and C. A. Meriles, “Detection and modeling of hole capture by single point defects under variable electric fields,” Nano Letters (2023), 10.1021/acs.nanolett.3c00860.
  * Manson and Harrison (2005) N. Manson and J. Harrison, “Photo-ionization of the nitrogen-vacancy center in diamond,” Diamond and Related Materials 14, 1705–1710 (2005).
  * Aslam  _et al._ (2013) N. Aslam, G. Waldherr, P. Neumann, F. Jelezko, and J. Wrachtrup, “Photo-induced ionization dynamics of the nitrogen vacancy defect in diamond investigated by single-shot charge state detection,” New Journal of Physics 15, 013064 (2013).
  * Brown  _et al._ (2019) K. J. Brown, E. Chartier, E. M. Sweet, D. A. Hopper, and L. C. Bassett, “Cleaning diamond surfaces using boiling acid treatment in a standard laboratory chemical hood,” Journal of Chemical Health & Safety 26, 40–44 (2019).
  * Falk  _et al._ (2013) A. L. Falk, B. B. Buckley, G. Calusine, W. F. Koehl, V. V. Dobrovitski, A. Politi, C. A. Zorman, P. X.-L. Feng, and D. D. Awschalom, “Polytype control of spin qubits in silicon carbide,” Nature Communications 4, 1819 (2013).
  * Diler  _et al._ (2020) B. Diler, S. J. Whiteley, C. P. Anderson, G. Wolfowicz, M. E. Wesson, E. S. Bielejec, F. J. Heremans, and D. D. Awschalom, “Coherent control and high-fidelity readout of chromium ions in commercial silicon carbide,” npj Quantum Information 6 (2020), 10.1038/s41534-020-0247-7.
  * Zhang  _et al._ (2020) G. Zhang, Y. Cheng, J. P. Chou, and A. Gali, “Material platforms for defect qubits and single-photon emitters,” Applied Physics Reviews 7 (2020), 10.1063/5.0006075.
  * Spindlberger  _et al._ (2019) L. Spindlberger, A. Csóré, G. Thiering, S. Putz, R. Karhu, J. U. Hassan, N. T. Son, T. Fromherz, A. Gali, and M. Trupke, “Optical properties of vanadium in 4 h silicon carbide for quantum technology,” Physical Review Applied 12 (2019), 10.1103/PhysRevApplied.12.014015.


# A two-dimensional programmable tweezer array of fermions

Zoe Z. Yan1,∗, Benjamin M. Spar1,∗, Max L. Prichard1, Sungjae Chi1, Hao-Tian
Wei2,3, Eduardo Ibarra-García-Padilla2,3, Kaden R. A. Hazzard2,3, Waseem S.
Bakr1 1 Department of Physics, Princeton University, Princeton, New Jersey
08544, USA
2 Department of Physics and Astronomy, Rice University, Houston, Texas 77005,
USA
3 Rice Center for Quantum Materials, Rice University, Houston, Texas 77005,
USA

###### Abstract

We prepare high-filling two-component arrays of up to fifty fermionic atoms in
optical tweezers, with the atoms in the ground motional state of each tweezer.
Using a stroboscopic technique, we configure the arrays in various two-
dimensional geometries with negligible Floquet heating. Full spin- and
density-resolved readout of individual sites allows us to post-select near-
zero entropy initial states for fermionic quantum simulation. We prepare a
correlated state in a two-by-two tunnel-coupled Hubbard plaquette,
demonstrating all the building blocks for realizing a programmable fermionic
quantum simulator.

††preprint: APS/123-QED Figure 1:  2D stroboscopic tweezer technique and
lifetimes. (a) Two crossed acousto-optic modulators spaced in a 4$f$
configuration generate the array. (b) Illustration of the principle of
stroboscopic array generation of an 8-site ring. For a strobe frequency
$f_{s}$, each column of the array is turned on for a quarter of the period
$1/f_{s}$, generating a time-averaged potential shown in (c). (d) Lifetime of
an atom in the ground vibrational state of a tweezer versus strobe frequency,
with the red point at 0 kHz indicating the non-strobed lifetime. The dashed
line shows the theoretical prediction, and grey shading indicates the
systematic uncertainties on the tweezer waist. The inset shows an example of a
decay curve of population in the ground state for $f_{s}\,{=}\,513$ kHz with
an exponential fit.

Ultracold atoms in optical tweezer arrays have become a popular platform for
quantum simulation, computation, and metrology [1]. The tweezer platform has
recently witnessed rapid breakthroughs, ranging from the development of
precise optical clocks [2, 3] to the demonstration of entangling operations
[4, 5, 6, 7]. The realization of defect-free arbitrary geometries [8, 9], in
particular in two dimensions, has paved the way for studying rich quantum
many-body physics with localized Rydberg atoms, including frustrated spin
models on a triangular lattice [10, 11], topological phases in a zig-zag chain
[12], and quantum spin liquids with atoms placed on the links of a kagome
lattice [13].

The versatility of tweezer arrays has also been extended to systems of
itinerant atoms where quantum statistics play a role [14, 15, 16, 17, 18, 19].
In particular, tunnel-coupled arrays have been realized for small systems of
bosonic [14] and fermionic [15, 16, 17, 18] atoms in one dimensional arrays.
If these experiments can be scaled, they would constitute a bottom-up approach
toward quantum simulation that complements optical lattice experiments with
quantum gas microscopes, which currently lie at the forefront of studying one-
and two-dimensional Hubbard models [20, 21, 22, 23, 24, 25, 26, 27, 28, 29].
The difficulty of reconfiguring microscope experiments has led to an almost
exclusive focus on physics in square lattices (Ref. [30] is a recent
exception). Programmable Hubbard tweezer arrays would allow the extension of
site-resolved studies to arbitrary lattice geometries that bring additional
ingredients into play, including frustration, topology, and flat-band physics.

Hubbard tweezer arrays also have the potential to address another major
challenge for optical lattice experiments: the preparation of low-entropy
phases of fermions. In optical lattice experiments, the entropy of the gas is
limited by evaporative cooling, which is hindered by poor efficiencies at low
temperatures. Entropy redistribution schemes relying on the flow of entropy
away from gapped phases have been proposed [31, 32] and experimentally
explored [33], but they have not resulted in significant reduction of achieved
temperatures.

Here we show that stroboscopic optical tweezer arrays can be used to prepare
low entropy fermionic systems with arbitrary two-dimensional (2D) geometry.
The low entropy is possible due to several features particular to this
platform. First, in loading a tweezer from a degenerate Fermi gas, the tweezer
acts as a “dimple trap,” wherein the local Fermi temperature
($T_{\mathrm{F}}$) is significantly higher than in the bulk gas. Since the
fraction of atoms loaded into the tweezers is low, the temperature of the
system remains approximately fixed to the bulk gas’s temperature, but the
tweezers’ phase space density is enhanced. Furthermore, the occupancy of the
lowest level of each tweezer (given by the Fermi-Dirac distribution) is close
to unity. This enables the preparation of a state with two atoms in the ground
motional state (one per spin state) on every tweezer with high fidelity, as
first demonstrated in Ref. [34]. Second, the system can be evolved from the
band insulator into a correlated state via an adiabatic ramp-on of additional
sites, taking advantage of independent tunability of each lattice site. We
have previously shown that this technique can be used to prepare a state with
antiferromagnetic correlations in an eight-site Fermi Hubbard chain [18]. We
extend this approach to 2D and show that any pre-ramp entropy in the system
can be effectively eliminated by post-selection on the atom number in each
spin state. Post-selection is enabled by spin- and density-resolved readout
[22, 35], which we implement in a bilayer imaging scheme.

Figure 2:  Examples of band insulators of different geometries, showing (a)
rectangular $5\times 5$, (2) 21-site Lieb plaquette, (3) $4\times 5$
triangular, and (4) octagonal ring arrays. Only $|{\uparrow}\rangle$atoms are
imaged and the sites here are not tunnel-coupled. The top row shows single
shots with perfect filling of the $|{\uparrow}\rangle$ state, and the bottom
row shows average images. Deviations of the atom positions in the single-shot
images are due to quantization onto the lattice for imaging. Average fillings
of $|{\uparrow}\rangle$ are (93, 92, 91, 89)%, accounting for imaging fidelity
of 98.5%, out of (411, 254 , 275, 100) shots.

The experimental cycle, including tweezer loading, is the same as detailed in
Ref. [18]. Tweezers are loaded from a bulk Fermi gas at $T/T_{\rm
F}\,{\approx}\,0.2$ that is a balanced mixture of the lowest and third lowest
hyperfine ground states ($|{\uparrow}\rangle$, $|{\downarrow}\rangle$,
respectively). Our scheme for generating 2D arrays uses two crossed acousto-
optical modulators [Fig. 1(a)]. The tweezers are produced using light with a
wavelength of $780$ nm, and their waist at the atoms is $1000^{+180}_{-140}$
nm. Radiofrequency tones for both AOMs are generated by a two-channel
arbitrary waveform generator, with a tone separation of 8 MHz corresponding to
a tweezer spacing of 1350 nm in the atom plane. The aperture size and
bandwidth of the modulators currently limit us to ${\sim}\,9$ tweezers in each
direction. The beat frequency of neighboring tweezers is $>100$ times larger
than typical tweezer depths, leading to negligible parametric heating.

Homogenizing the tweezer depths is particularly challenging for 2D arrays
generated using crossed AOMs. A common approach used in Rydberg tweezer
experiments is to apply a static set of frequencies consisting of $n_{x}$ and
$n_{y}$ tones for the $x$\- and $y$\- directional AOMs, respectively. This
generates a rectangular array of $n_{x}n_{y}$ sites; however, the
$n_{x}\,{+}\,n_{y}$ degrees of freedom from the signal strength of each tone
are insufficient to independently tune the depth of each tweezer. Better
homogeneity can be achieved by tuning the relative phases of the tones, but
the typical resultant inhomogeneity still exceeds 1%. Tunnel-coupled arrays
have more stringent requirements for homogeneity, since the energy offsets in
tweezers of typical depth $\sim h{\times}50\,$kHz must be controlled to within
tunneling energies of $\sim h{\times}250\,$Hz, or better than 0.5%.

To homogenize arrays to within this precision and produce arrays with
arbitrary geometry, we introduce a stroboscopic tweezer technique. We generate
the array one column at a time, with different $y-$directional tones applied
in every timestep [Fig. 1(b)]. Effectively, the atoms experience a time-
averaged potential of concatenated 1D arrays, as long as the strobe rate
$f_{s}$ far exceeds the tweezers’ harmonic trap frequencies. As the typical
axial (radial) trap frequencies are around 2.5 (15) kHz, we need strobe rates
over an order of magnitude higher to avoid significant Floquet heating of the
atoms.

We verify that the stroboscopic scheme is compatible with long lifetimes in
the tweezer ground vibrational state with the following study. We measure the
dependence of the lifetime in the lowest vibrational state on $f_{s}$ in a
two-site strobed array, varying the strobe rate from 163-1083 kHz [36]. Higher
frequencies are inaccessible due to limitations on the AOM response rate, set
by the speed of sound and beam size in the crystal. We also compare the
lifetimes to that of a static (non-strobed) tweezer, which is limited by
background gas collisions and off-resonant photon scattering due to the
trapping light. Consistent with expectations, the lowest strobe rates give the
shortest lifetimes in the ground state [Fig. 1(d)]. Measurements and numerics
using a discrete variable representation (DVR) method [37, 38, 36] both
indicate that Floquet heating decreases exponentially with increasing $f_{s}$
and is negligible for $f_{s}\gtrsim 250$ kHz, although the numerics
underestimate the threshold frequency range below which severe heating occurs
by ${\sim}\,18\%$.

We demonstrate loading the arrays with band insulators of fermions with high
fidelity using the stroboscopic method. These band insulators serve as low
entropy initial states for fermionic quantum simulation. As proofs-of-
principle, we implement a rectangular $5\times 5$ array, 21-site Lieb
plaquette, triangular $4\times 5$ array, and an 8-site octagonal ring (Fig. 2)
[36]. The tweezers are homogenized using a density balancing algorithm where
the number of required experimental shots is almost independent of the array
size [18]. In these examples, the sites are not tunnel-coupled due to the
large separations. Readout is accomplished by transferring the atoms into a 2D
lattice of 752 nm spacing, which oversamples the tweezer array, and performing
Raman sideband cooling on the $|{\uparrow}\rangle$ atoms after removal of the
$|{\downarrow}\rangle$ atoms [23, 18] with a detection fidelity of 98.5%.
Throughout these different geometries, the loading fidelity of a single spin
averages to 92%/site, corrected for imaging infidelity, indicating a low
entropy of loading in the array. As in previous work [34, 18], the tweezer
depths are chosen so the predominant type of defect in each tweezer is a
missing particle rather than an extra one in a higher motional state. In these
data, we only measure one of the spin states in a given experimental shot, due
to the problem of light-assisted collisions, which necessitates the removal of
the other spin state before imaging [39].

Figure 3:  Bilayer imaging procedure and entropy reduction through post-
selection. (a) Atoms in $|{\uparrow}\rangle$ (yellow) and
$|{\downarrow}\rangle$ (blue) are initially trapped in the tweezers, then
adiabatically loaded into (b) a 2D lattice with vertical waist of 75 $\mu$m,
where $|{\uparrow}\rangle$ is transferred to $|{\tilde{\uparrow}}\rangle$
(red). (c) A magnetic field gradient is applied to separate the spins in the
vertical direction, after which (d) two lightsheet potentials turn on to fix
the $z$-positions. (e) The lightsheets are further separated to 25 $\mu$m
separation. Raman sideband imaging commences, producing simultaneous images of
both spin states. (f) shows a single shot image of $|{\downarrow}\rangle$ and
$|{\tilde{\uparrow}}\rangle$ originally from a $3\times 5$ rectangular array.
(g) Probability distribution versus number of atoms in each spin state over
972 shots. Here, all images with doublons (65 shots) were not used. (h) By
post-selecting on the maximum number of holes, effective entropy can be
reduced by varying amounts.

To circumvent this problem and obtain full density- and spin-resolution, we
adopt a high-fidelity bilayer imaging scheme [40, 35, 41, 42], which also
allows the reduction of entropy upon post-selection. Bilayer density- and
spin-readout was first accomplished in fermionic quantum gas microscope
experiments in a superlattice charge-pumping scheme [35]. Our method is
conceptually similar but involves no superlattice (Fig. 3). Starting with
tweezer-trapped atoms, we adiabatically turn off the tweezer and turn on a 2D
lattice of 1064 nm to 60 $E_{\mathrm{R}}$ and a vertical trap frequency of 1.2
kHz in 5 ms. The magnetic field is brought to 572 G, where we perform a spin-
flip of $|{\uparrow}\rangle$ to the second-lowest hyperfine state,
$|{\tilde{\uparrow}}\rangle$, with an efficiency exceeding 99%, and then
decrease the field to near 0 G. Atoms in $|{\tilde{\uparrow}}\rangle$ and
$|{\downarrow}\rangle$ have a greater differential magnetic moment than those
in $|{\uparrow}\rangle$, enabling the Stern-Gerlach separation of these
populations to ${\sim}\,9$$\mu$m using a $z$-magnetic gradient of 168 G/cm in
the 2D lattice at a depth of 280 $E_{\mathrm{R}}$. We turn on two lightsheet
potentials [23]–highly anisotropic beams, each with $z$-directional trap
frequencies of 26 kHz–and linearly ramp their vertical separation to 25 $\mu$m
for imaging. We measure a combined transport and spin identification fidelity
of 98.7%. Finally, we image the atoms using Raman sideband cooling
simultaneously in both layers, with the 2D lattice depth at 2500
$E_{\mathrm{R}}$ and the two lightsheet $z$-trap frequencies at 70 kHz.
Resulting fluorescence is collected by a high numerical aperture objective
with atoms in the two planes focused onto two different active areas of a CCD
camera. Imaging fidelity is 98% (97%) for the layer of
$|{\tilde{\uparrow}}\rangle$ ($|{\downarrow}\rangle$) atoms.

Bilayer imaging enables reduction of the effective entropy associated with the
initial state of the tweezer array (the band insulator) through post-
selection. The initial entropy per particle of the tweezer ensemble, assuming
independent tweezers and single-band occupation, is given by

$\frac{S}{\langle N\rangle}=-\frac{k_{\mathrm{B}}}{p}\big{(}p\log
p+(1-p)\log(1-p)\big{)}.$ (1)

where $p$ is the probability to load one spin on a site. With a typical
loading efficiency of $p\,{=}\,0.907(3)$, the array starts with 0.34(1)
$k_{\mathrm{B}}$ per particle, with entropy entering from microstates with
undesired holes. By selecting only images with the population per spin state
equal to the number of loading tweezers, we can effectively choose a subsample
with $S\,{=}\,0$. Importantly, this post-selection criterion eliminates the
initial state entropy even after changing the filling of the system (by
introducing additional tweezers) to prepare a correlated state. The post-
selection criterion can be relaxed to use more images from the experiment at
the cost of introducing additional initial state entropy. This tradeoff is
illustrated in Figs. 3(g-h) for a $3\times 5$ array in which
$|{\uparrow}\rangle$ and $|{\downarrow}\rangle$ had average $p\,{=}\,$0.914(3)
and 0.900(3), respectively (not accounting for imaging fidelity). Out of 972
images, 12% had perfect filling of 15 fermions of each species. However, even
keeping images with up to two holes, or over 50% of shots, still results in a
low entropy of 0.17(1) $k_{\mathrm{B}}$ per particle, which is favorable
compared to state-of-the-art optical lattice experiments that range from
0.3-0.5 $k_{\mathrm{B}}$ per particle [43, 25, 44].

While post-selection can be used to reduce the effective entropy of the
initial state to near zero, subsequent ramps to correlated states will
inevitably introduce additional entropy. Numerical simulations of the
dynamical ramps in small systems indicate this extra entropy should be low for
defect-free initial configurations. For example, for the ramp used in our
previous work with an eight-site chain [18], the ramp is expected to have
introduced an additional entropy of 0.04 $k_{\mathrm{B}}$ per particle when
starting with a defect-free state, but the presence of even a single localized
hole would lead to a significant entropy increase of 0.2-0.3 $k_{\mathrm{B}}$
per particle depending on the position of the hole. The entropy reported in
Fig. 3(h) should therefore be treated only as a lower bound for future
experiments.

Post-selection on spin and density in this context should be distinguished
from the context of optical lattice-based quantum gas microscopy measurements.
For example, in a recent study with a fermionic microscopes [44], spin- and
density- readout enabled post-selection of half-filled systems with zero total
magnetization, keeping ${\sim}\,$9% of data. However, post-selection there did
not eliminate the finite spin temperature in the initial state.

Figure 4:  Low entropy preparation of a $2\times 2$ array. a) We load two
atoms per site on one diagonal of the array. b) We create a correlated state
by ramping on the additional lattice sites and increasing the scattering
length to introduce on-site interactions. For the following data, we work with
$t_{x}=140(5)$ Hz, $t_{y}=220(5)$ Hz and $U/\bar{t}=3.4(2)$. c) Measured spin-
spin correlations enabled by the bilayer imaging scheme. d) Best fit (top) and
measured (bottom) microstate populations for 671 post-selected experimental
shots. The fit gives an entropy in the range [0,0.09] $k_{B}$ per particle.
Insets are shown for the two most common states.

Equipped with the ability to load near-zero-entropy band insulators after
post-selection, we implement the simplest building block of a two-dimensional
Fermi-Hubbard model: a tunnel-coupled $2\times 2$ plaquette. The single-band
Hamiltonian is

$\displaystyle\hat{H}$ $\displaystyle=-\sum_{\langle
i,j\rangle_{x},\sigma}t_{x}(\hat{c}^{{\dagger}}_{i\sigma}\hat{c}_{j\sigma}{+}{\rm
h.c.})\,{-}\sum_{\langle
i,j\rangle_{y},\sigma}t_{y}(\hat{c}^{{\dagger}}_{i\sigma}\hat{c}_{j\sigma}{+}{\rm
h.c.})$
$\displaystyle\,{+}\,\sum_{i}U_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}{+}\sum_{i,\sigma}\Delta_{i}\hat{n}_{i\sigma},$
(2)

where $\hat{c}^{{\dagger}}_{i\sigma}$ is the fermionic creation operator of
spin $\sigma$ at site $i$, $\hat{n}_{i\sigma}$ is the number operator,
$t_{x(y)}$ is the tunneling matrix element in the $x$($y$) direction, $\Delta$
is the energy offset, and $U$ is on-site interaction between opposite spin
states. We start by loading two diagonal sites in a rectangular array with
vertical (horizontal) spacing of 1520 (1690) nm [Fig. 4 (a,b)]. The correlated
state at half-filling is prepared by adiabatically ramping on the two opposing
diagonal sites in 50 ms [18], with tunnelings of
$t_{x}\,[t_{y}]=h{\times}140(5)[220(5)]$ Hz in the final configuration [36].
We also ramp $U/\bar{t}$ from 0 to 3.4(2) in the same time using the Feshbach
resonance. Here, $\bar{t}\,{=}\,(t_{x}+t_{y})/2$.

The resulting spin-spin correlations are shown in Fig. 4(c), which depicts
$C_{ij}\,{=}\langle S_{z,i}S_{z,j}\rangle{-}\langle S_{z,i}\rangle\langle
S_{z,j}\rangle$, where
$S_{z,i}{\equiv}\frac{1}{2}(n_{\uparrow,i}{-}n_{\downarrow,i})$. Here, data
were post-selected to include only images that contained two
$|{\uparrow}\rangle$ and two $|{\downarrow}\rangle$ atoms, for a total of 673
experimental cycles. With full spin- and density-readout, we are able
reconstruct the diagonals of the density matrix
$\rho\,{=}\,|\Psi\rangle\langle\Psi|$ in the basis of allowed number states
(with Hilbert space size ${{4}\choose{2}}^{2}{=}\,36$), and compare data with
theory. In Fig. 4(d), we plot the experimental population in each microstate
together with the populations expected theoretically for the plaquette ground
state, which are consistent within error bars. Here, we reduce our statistical
errors by taking advantage of the spin-symmetry of the Hubbard Hamiltonian to
average the probabilities for spin-reversed microstates. Furthermore, we fit
the temperature of the canonical ensemble to best reproduce the distribution
of microstates. The fit gives an upper bound for the temperature of
$k_{B}T\sim 0.3\bar{t}$ (with the fit losing sensitivity below that
temperature). This corresponds to an entropy in the range [0,0.09] $k_{B}$ per
particle, which is consistent with the prediction from simulating the ramp
dynamics (entropy gain of 0.02 $k_{B}$ per particle).

In conclusion, we have realized a 2D tweezer array of fermions with software-
programmable geometry using a novel stroboscopic technique that allows
independent control over all tweezer depths and positions. We have realized
the building blocks to implement programmable 2D Fermi-Hubbard models, and
demonstrated these on a small scale. Future work will focus on increasing the
system size of the tunnel-coupled arrays. A natural target for future work
will be few-leg ladder systems. For example, two-leg triangular ladder systems
can be used to explore the $J_{1}-J_{2}$ model, including the special case of
the Majumdar-Ghosh model and its valence-bond solid ground states [45].
Furthermore, upon introducing spin-imbalance and hole-doping, a triangular
two-leg ladder is predicted to host magnon-hole binding at energy scales set
by the tunneling, rather than the superexchange [46]. Multi-leg triangular
ladders may potentially host other exotic states such as a chiral spin liquid
at half filling and intermediate $U/t$ that evolves to a
120∘-antiferromagnetic order at strong $U/t$ [47]. Ultimately, fully 2D
tunnel-coupled arrays with arbitrary geometry will be a rich playground for
exploring novel phases of correlated fermions.

We would like to thank Grace Sommers and Elmer Guardado-Sanchez for
experimental assistance and Martin Zwierlein’s group for the loan of
equipment. The experimental work was supported by the NSF (grant no. 2110475),
the David and Lucile Packard Foundation (grant no. 2016-65128), and the ONR
(grant no. N00014-21-1-2646). K.R.A.H, H.-T.W., and E.I.-G.-P. acknowledge
support from the Welch Foundation through Grant No. C1872, the Office of Naval
Research Grant No. N00014-20-1-2695, and the National Science Foundation
through Grant No. PHY1848304. K.R.A.H. also benefited from discussions at the
KITP, which was supported in part by the National Science Foundation under
Grant No. NSF PHY-1748958.

∗ These authors contributed equally to this work.

## References

  * Kaufman and Ni [2021] A. M. Kaufman and K.-K. Ni, Nat. Phys. 17, 1324 (2021).
  * Madjarov  _et al._ [2019] I. S. Madjarov, A. Cooper, A. L. Shaw, J. P. Covey, V. Schkolnik, T. H. Yoon, J. R. Williams, and M. Endres, Phys. Rev. X 9, 041052 (2019).
  * Young  _et al._ [2020] A. W. Young, W. J. Eckner, W. R. Milner, D. Kedar, M. A. Norcia, E. Oelker, N. Schine, J. Ye, and A. M. Kaufman, Nature (London) 588, 408 (2020).
  * Graham  _et al._ [2019] T. M. Graham, M. Kwon, B. Grinkemeyer, Z. Marra, X. Jiang, M. T. Lichtman, Y. Sun, M. Ebert, and M. Saffman, Phys. Rev. Lett. 123, 230501 (2019).
  * Levine  _et al._ [2019] H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, and M. D. Lukin, Phys. Rev. Lett. 123, 170503 (2019).
  * Omran  _et al._ [2019] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, _et al._ , Science 365, 570 (2019).
  * Madjarov  _et al._ [2020] I. S. Madjarov, J. P. Covey, A. L. Shaw, J. Choi, A. Kale, A. Cooper, H. Pichler, V. Schkolnik, J. R. Williams, and M. Endres, Nat. Phys. 16, 857 (2020).
  * Endres  _et al._ [2016] M. Endres, H. Bernien, A. Keesling, H. Levine, E. R. Anschuetz, A. Krajenbrink, C. Senko, V. Vuletic, M. Greiner, and M. D. Lukin, Science 354, 1024 (2016).
  * Barredo  _et al._ [2016] D. Barredo, S. De Léséleuc, V. Lienhard, T. Lahaye, and A. Browaeys, Science 354, 1021 (2016).
  * Lienhard  _et al._ [2018] V. Lienhard, S. de Léséleuc, D. Barredo, T. Lahaye, A. Browaeys, M. Schuler, L.-P. Henry, and A. M. Läuchli, Phys. Rev. X 8, 021070 (2018).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, _et al._ , Nature 595, 233 (2021).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Science 365, 775 (2019).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, _et al._ , Science 374, 1242 (2021).
  * Kaufman  _et al._ [2014] A. M. Kaufman, B. J. Lester, C. M. Reynolds, M. L. Wall, M. Foss-Feig, K. R. A. Hazzard, A. M. Rey, and C. A. Regal, Science 345, 306 (2014).
  * Murmann  _et al._ [2015] S. Murmann, A. Bergschneider, V. M. Klinkhamer, G. Zürn, T. Lompe, and S. Jochim, Phys. Rev. Lett. 114, 080402 (2015).
  * Bergschneider  _et al._ [2019] A. Bergschneider, V. M. Klinkhamer, J. H. Becher, R. Klemt, L. Palm, G. Zürn, S. Jochim, and P. M. Preiss, Nature Physics 15, 640 (2019).
  * Becher  _et al._ [2020] J. H. Becher, E. Sindici, R. Klemt, S. Jochim, A. J. Daley, and P. M. Preiss, Phys. Rev. Lett. 125, 180402 (2020).
  * Spar  _et al._ [2021] B. M. Spar, E. Guardado-Sanchez, S. Chi, Z. Z. Yan, and W. S. Bakr, arXiv preprint arXiv:2110.15398 (2021).
  * Young  _et al._ [2022] A. W. Young, W. J. Eckner, N. Schine, A. M. Childs, and A. M. Kaufman, arXiv:2202.01204 (2022).
  * Cheuk  _et al._ [2016] L. W. Cheuk, M. A. Nichols, K. R. Lawrence, M. Okan, H. Zhang, E. Khatami, N. Trivedi, T. Paiva, M. Rigol, and M. W. Zwierlein, Science 353, 1260 (2016).
  * Parsons  _et al._ [2016] M. F. Parsons, A. Mazurenko, C. S. Chiu, G. Ji, D. Greif, and M. Greiner, Science 353, 1253 (2016).
  * Boll  _et al._ [2016] M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Science 353, 1257 (2016).
  * Brown  _et al._ [2017] P. T. Brown, D. Mitra, E. Guardado-Sanchez, P. Schauß, S. S. Kondov, E. Khatami, T. Paiva, N. Trivedi, D. A. Huse, and W. S. Bakr, Science 357, 1385 (2017).
  * Nichols  _et al._ [2019] M. A. Nichols, L. W. Cheuk, M. Okan, T. R. Hartke, E. Mendez, T. Senthil, E. Khatami, H. Zhang, and M. W. Zwierlein, Science 363, 383 (2019).
  * Brown  _et al._ [2019] P. T. Brown, D. Mitra, E. Guardado-Sanchez, R. Nourafkan, A. Reymbaut, C.-D. Hébert, S. Bergeron, A.-M. Tremblay, J. Kokalj, D. A. Huse, _et al._ , Science 363, 379 (2019).
  * Vijayan  _et al._ [2020] J. Vijayan, P. Sompet, G. Salomon, J. Koepsell, S. Hirthe, A. Bohrdt, F. Grusdt, I. Bloch, and C. Gross, Science 367, 186 (2020).
  * Ji  _et al._ [2021] G. Ji, M. Xu, L. H. Kendrick, C. S. Chiu, J. C. Brüggenjürgen, D. Greif, A. Bohrdt, F. Grusdt, E. Demler, M. Lebrat, and M. Greiner, Phys. Rev. X 11, 021022 (2021).
  * Koepsell  _et al._ [2019] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Nature (London) 572, 358 (2019).
  * Bohrdt  _et al._ [2021] A. Bohrdt, L. Homeier, C. Reinmoser, E. Demler, and F. Grusdt, Annals of Physics 435, 168651 (2021), special issue on Philip W. Anderson.
  * Yang  _et al._ [2021] J. Yang, L. Liu, J. Mongkolkiattichai, and P. Schauss, PRX Quantum 2, 020344 (2021).
  * Bernier  _et al._ [2009] J.-S. Bernier, C. Kollath, A. Georges, L. De Leo, F. Gerbier, C. Salomon, and M. Köhl, Phys. Rev. A 79, 061601 (2009).
  * Lubasch  _et al._ [2011] M. Lubasch, V. Murg, U. Schneider, J. I. Cirac, and M.-C. Bañuls, Phys. Rev. Lett. 107, 165301 (2011).
  * Chiu  _et al._ [2018] C. S. Chiu, G. Ji, A. Mazurenko, D. Greif, and M. Greiner, Phys. Rev. Lett. 120, 243201 (2018).
  * Serwane  _et al._ [2011] F. Serwane, G. Zürn, T. Lompe, T. Ottenstein, A. Wenz, and S. Jochim, Science 332, 336 (2011).
  * Koepsell  _et al._ [2020] J. Koepsell, S. Hirthe, D. Bourgund, P. Sompet, J. Vijayan, G. Salomon, C. Gross, and I. Bloch, Phys. Rev. Lett. 125, 010403 (2020).
  * [36] See supplement.
  * Light and Carrington Jr [2000] J. C. Light and T. Carrington Jr, Advances in Chemical Physics 114, 263 (2000).
  * Wall  _et al._ [2015] M. L. Wall, K. R. A. Hazzard, and A. M. Rey, Phys. Rev. A 92, 013610 (2015).
  * Gross and Bakr [2021] C. Gross and W. S. Bakr, Nature Physics 17, 1316 (2021).
  * Preiss  _et al._ [2015] P. M. Preiss, R. Ma, M. E. Tai, J. Simon, and M. Greiner, Phys. Rev. A 91, 041602 (2015).
  * Hartke  _et al._ [2020] T. Hartke, B. Oreg, N. Jia, and M. Zwierlein, Phys. Rev. Lett. 125, 113601 (2020).
  * Gall  _et al._ [2021] M. Gall, N. Wurz, J. Samland, C. F. Chan, and M. Köhl, Nature 589, 40 (2021).
  * Mazurenko  _et al._ [2017] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature (London) 545, 462 (2017).
  * Sompet  _et al._ [2021] P. Sompet, S. Hirthe, D. Bourgund, T. Chalopin, J. Bibo, J. Koepsell, P. Bojović, R. Verresen, F. Pollmann, G. Salomon, C. Gross, T. Hilker, and I. Bloch, arXiv:2103.10421 (2021).
  * Girvin and Yang [2019] S. M. Girvin and K. Yang,  _Modern Condensed Matter Physics_ (Cambridge University Press, 2019).
  * Morera  _et al._ [2021] I. Morera, A. Bohrdt, W. W. Ho, and E. Demler, arXiv:2106.09600 (2021).
  * Szasz  _et al._ [2020] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Phys. Rev. X 10, 021042 (2020).
  * Neuhasuer and Baer [1989] D. Neuhasuer and M. Baer, The Journal of Chemical Physics 90, 4351 (1989), https://doi.org/10.1063/1.456646 .
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, SciPost. Phys. 7, 97 (2019).

Supplemental Information:
A two-dimensional programmable tweezer array of fermions

##  1 Measurement of Floquet heating

To measure the ground state lifetimes shown in Fig. 1 of the main text, we
perform the following experiment. We load two tweezers, spaced 8.1 $\mu$m
apart to avoid any overlap, each with a spin up and a spin down atom. Atoms in
higher vibrational states are removed by lowering the tweezer depth and
applying a magnetic field gradient, as described in [18]. The tweezer depth is
increased to 50 kHz, and a variable hold time is applied. The tweezer light is
strobed such that only one tweezer is on at a time, resulting in Floquet
heating for low strobe frequencies. Finally, atoms in newly populated higher
vibrational states are again removed with an identical spilling process, and
the remaining $|{\uparrow}\rangle$ atoms are imaged.

##  2 Theoretical calculations of Floquet heating

We determine the Floquet heating from the strobing of the trap by directly
calculating the dynamics of a single atom in a tweezer potential using a
method based on a discrete variable representation (DVR) of Hilbert space
[37], which was applied to ground state properties of tunnel-coupled tweezers
in [38, 14]. The tweezer potential is

$V(\vec{r})=-\frac{V_{0}}{1+\frac{z^{2}}{z_{R}^{2}}}\exp\left[-\frac{2r^{2}}{w_{0}^{2}(1+\frac{z^{2}}{z_{R}^{2}})}\right]$
(S1)

where $V_{0}$ is the trap depth when the tweezer is on, $z_{R}$ is the
Rayleigh range, and $w_{0}$ is the trap waist, and we use the measured
parameters. We use a DVR basis of spatial-parity-adapted sinc functions on a
three-dimensional cubic grid of $(N_{x},N_{y},N_{z})$ points spanning
$(0,0,0)$ to $(L_{x},L_{y},L_{z})$ to represent our wavefunctions. Only even
parity states occur for the dynamics considered.

The state at time $t$ after the strobed time evolution is

$\displaystyle\ket{\psi(t)}=\underbrace{\left[e^{-i{\hat{T}}\delta
t}e^{-i({\hat{T}}+{\hat{V}})\delta t}\right]\cdots\left[e^{-i{\hat{T}}\delta
t}e^{-i({\hat{T}}+{\hat{V}})\delta t}\right]}_{t/2\delta t\text{
times}}\ket{i}$ (S2)

where ${\hat{T}}=-\hbar^{2}\nabla^{2}/2m+V_{\text{abs}}(\vec{r})$ is the
kinetic energy operator ($\hbar$ is Planck’s constant, and $m$ is the mass of
6Li) plus a potential ${\hat{V}}_{\text{abs}}$ to handle boundary conditions
discussed below, ${\hat{V}}$ is the potential energy operator associated with
Eq. (S1), and $\delta t=1/(2f_{s})$, and assuming $t/2\delta t$ is an integer
(the micromotion between such times is discussed later). The initial state
$\ket{i}$ is assumed to be the ground state of the system in the time-averaged
potential, i.e. of $H_{1/2}=T+V/2$. Directly solving Eq. (S2) is challenging
due to three factors: (1) the tweezer potential is 3D, (2) it involves a wide
range of length scales (the wavefunction localization length, the Gaussian
waist, and the Rayleigh range), and (3) the dynamics problem must account for
a huge separation of timescales when $f_{s}$ is large, the short time of the
pulses and the long lifetime of the atoms, a ratio of timescales of more than
$10^{7}$.

We solve this by rewriting Eq. (S2). First we can diagonalize to obtain
$e^{-i{\hat{T}}\delta t}=U_{0}D_{0}U_{0}^{-1}$ where $D_{0}$ is diagonal and
$U_{0}$ is the matrix whose columns are eigenvectors of ${\hat{T}}$, obtained
by first diagonalizing $\hat{T}$ and then exponentiating. ($U_{0}$ need not be
unitary since ${\hat{T}}$ has non-Hermitian terms for the absorbing
potential.) Similarly $e^{-i({\hat{T}}+{\hat{V}})\delta
t}=U_{1}D_{1}U_{1}^{-1},$. Therefore the terms in brackets – denote it
$M\equiv e^{-i{\hat{T}}\delta t}e^{-i({\hat{T}}+{\hat{V}})\delta t}$ – in Eq.
(S2) can be rewritten

$M=U_{0}D_{0}U_{0}^{-1}U_{1}D_{1}U_{1}^{-1}.$ (S3)

One can multiply these matrices, and diagonalize the result to find
$M=U_{2}D_{2}U_{2}^{-1}$, so

$\ket{\psi(t)}=M^{t/2\delta t}\ket{i}=U_{2}D_{2}^{t/2\delta
t}U_{2}^{-1}\ket{i}.$ (S4)

The micromotion at times between integer multiples of $2\delta t$ can be
treated by propagating to the largest multiple of $1/f_{s}$ as above, and then
calculating the time evolution within a single period, which is
straightforward and efficient in terms of the diagonalized evolution matrices
($D_{a}$, $U_{a}$) already obtained. The computational cost is dominated by
diagonalizing to find $D_{0}$, $D_{1}$, and $D_{2}$, with several additional
multiplications. Although this involves full diagonalization of three large
matrices, it avoids the issues associated with separation of timescales.

To obtain accurate results, one must use a a fine enough DVR grid and a large
enough system. Here, a challenge arises that is absent for ground state
calculations: atoms can be excited from the ground state to scattering states
and thus escape from the trap. They may move rapidly, and thus would require
$L_{a}$ that grow linearly with time to capture, in practice many orders of
magnitude larger than we can calculate. However, in the experiment once
particles move sufficiently far outside the trap, the probability they return
is negligible. This can be accurately reproduced in the calculations by
including an absorbing boundary region in our simulations, a standard
technique in calculations of, e.g., chemical dynamics [48]. We use a potential

$V_{\text{abs}}(\vec{r})=-i\Gamma\\!\sum_{a\in
x,y,z}\frac{|a|-L_{a}^{(0)}}{L_{a}-L_{a}^{(0)}}.$ (S5)

This is zero inside the cubic region from the origin to
$(L_{x}^{(0)},L_{y}^{(0)},L_{z}^{(0)})$, and linearly increases as one moves
outside this region. The whole calculation is performed inside a box of size
$(L_{x},L_{y},L_{z})$ with $(N_{x},N_{y},N_{z})$ grid points. The calculations
in the main text are performed with
$(L_{x}^{(0)},L_{y}^{(0)},L_{z}^{(0)})=(3,3,7.2)w_{0}$, $N=(27,27,23)$;
$\Gamma=$57\text{\,}\mathrm{kHz}$$; and $L_{a}=L_{a}{(0)}+S_{a}$ where
$S_{a}\approx w_{0}$, and are well-converged in the number of grid points and
system size (Fig. S1). ($S_{a}$ varies slightly in the right panel: because we
fix the grid spacing and ensure that the last grid point in the non-absorbing
region doesn’t change with $L_{a}$, we choose $S_{a}$ closest to $w_{0}$ while
constraining $L_{a}$ to correspond to the location of the last grid points
along the $a$’th direction.) The $\Gamma$ and $S$ are chosen to be
sufficiently large so that particles leaving the trap are absorbed before
reaching the boundary of the calculation, but $\Gamma$ is maintained
sufficiently small so that there is no Zeno effect that would reflect the
particles. The choices we make are consistent with these conditions, as given
in Ref. [48], and our results are insensitive to the specific choice of the
parameters. For example, varying $\Gamma$ over a range of a factor of hundred
has a negligible effect on the atom lifetime. Calculated lifetimes are
converged to have visually negligible errors from the number of grid points,
system size, and absorbing boundary.

Using this method, we calculate the population of the ground state (of the
time-averaged potential) as a function of time. We find it gives a nearly
perfect exponential decay for the $f_{s}$ shown, with small deviations from
exponential at smaller $f_{s}$. The lifetimes shown in Fig. 1(d) of the main
text are obtained by fitting an exponential $e^{-t/\tau}$ to the obtained
dynamics. The numerical calculations do not include the effects that give the
finite lifetime in the unstrobed tweezer (background gas collisions and off-
resonant light scattering), and which are expected to be the reason for the
high-frequency lifetime saturation in Fig. 1(d). To incorporate those effects,
the theory curve in Fig. 1(d) is a simple interpolation
$1/\tau_{\text{eff}}=1/\tau+1/\tau_{\text{static}}$ where
$\tau_{\text{static}}$ is the lifetime in the unstrobed tweezer. The plot is
largely insensitive to the details of this interpolation.

Figure S1: Convergence of the calculated atom lifetime with respect to number
of DVR grid points (left) and system size (right). Left: Lifetimes for
different numbers of DVR grid points are converged well in the number of
gridpoints for $(N_{x},N_{y},N_{z})=(25,25,21)$ (or more) in systems of size
$(L^{(0)}_{x},L_{y}^{(0)},L_{z}^{(0)})=(3,3,7.2)w_{0}$. Right: lifetimes for
different system sizes are converged well for
$(L_{x},L_{y},L_{z})=(3,3,5.6)w_{0}$ and larger. Right panel is at fixed grid
point spacing $(0.16667,0.166667,0.4)w_{0}$, so the $N_{a}$ vary with $L_{a}$.
The spacing corresponds to $(N_{x},N_{y},N_{z})=(24,24,20)$ for the largest
system.

##  3 2D Array Generation

Our arrays are generated by two perpendicularly oriented AOMs from AA
Optoelectronic (MT110-B50A1,5-IR) driven by a two-channel arbitrary waveform
generator (Spectrum M4i.6621-x8 PCIe). The typical strobing scheme,
illustrated in Fig. 1 of the main text, generates one column of the array at
any given time. This gives the most flexibility for creating different
geometries. The two AOMs operate with the longitudinal acoustic mode with a
high speed of sound (4200 m/s) in the crystal. This speed, along with the
aperture size of $1.5\times 2$ mm2, sets the minimum possible dwell time on a
tweezer (the duration for which the tweezer is on). We find that dwell times
below 0.5 $\mu s$ produce tweezers with profiles that are significantly
distorted along the strobe axis. When generating a two-site array by strobing,
we measure that the tweezer waist along the strobe axis increases at 1.4% per
increase of 100 kHz strobing rate. Distortions in the tweezer profile can be
partly mitigated by applying a cosine-sum window function on the waveform that
is sent to the strobed AOM. Given the minimum dwell time per tweezer, the
strobe period $1/f_{s}$ grows linearly with the number of strobed columns.
This seemingly limits us to about 8 sites along the strobe axis due to Floquet
heating based on the results shown in Fig. 1 of the main text. However, we can
surpass this limitation with the technique discussed in the following
paragraph.

To reduce Floquet heating and distortion of the tweezer intensity profile, we
combine 2D multitone arrays (produced by driving both AOMs with a set of
radiofrequency tones) with strobing. For certain geometries, this allows us to
reduce tweezer dwell times by using smaller modulation depths. We define
modulation depth as the difference between the minimum and maximum relative
intensity used over one strobe period. The principle is illustrated in Fig. S2
for 50% modulation depth, and was applied for the rectangular, Lieb,
triangular, and ring lattices in Fig. 2 of the main text. For the example of
Fig. S2, a static 2D multitone configuration (modulation depth of 0%) can be
used in principle to create the geometry shown, but would not give enough
degrees of freedom to homogenize the depths of all sites.

The type of lattices that can be generated with this approach is restricted to
combinations of patterns that can be generated by using 2D multitone arrays.
Assume each of the two AOMs is driven by a set of tones such that it
individually produces tweezers at positions $P_{x},P_{y}$ respectively, where
$P_{x}\,{=}\,\\{p_{x,1},...p_{x,n}\\}$ and
$P_{y}\,{=}\,\\{p_{y,1},...p_{y,m}\\}$. The result when driving both AOMs will
be the Cartesian product of the two sets $P_{x}\times
P_{y}=\\{(p_{x,i},p_{y,j})\,|\,p_{x,i}\in P_{x},p_{y,j}\in P_{y}\\}$. For
example, triangular and Lieb lattices can be generated by switching between
two multitone rectangular arrays and overlaying a small amplitude strobed
array to homogenize the tweezer depths. In general, most 2D patterns of
interest that host some level of periodicity can be generated by switching
between two or three multitone rectangular arrays at full modulation depth and
overlaying a small amplitude strobed array for homogenization. When using 2D
multitone arrays, we avoid square geometries, as equidistant tweezer spacing
in $x$ and $y$ leads to diagonal sites having degenerate tones, which results
in low frequency beating between these tweezers and atom heating.

Figure S2:  Illustration of using 50% modulation depth to generate a $2\times
3$ array with strobe rate $f_{s}$, with (a) the individual strobe steps and
(b) the time-averaged potential. This procedure can be thought of as
superimposing a static $2\times 3$ grid on top of a 3-step array with 100%
modulation depth.

##  4 $2\times 2$ plaquette

We calibrate the tunnelings in the horizontal and vertical directions in the
$2\times 2$ array as follows. We first load atoms on only one site, and reduce
the tweezer intensities by 5 percent on all other sites except for the one in
the direction we want to measure the tunneling. At a non-interacting magnetic
field, we quickly initialize tunneling by zeroing the offsets among the target
sites [18], and we measure Rabi oscillations between the two sites [15, 14].

To calibrate the on-site interaction energy $U$, we perform radiofrequency
spectroscopy and measure the difference between single and double occupied
sites. However, the most sensitive way to measure $U/t$ is to fit the
distribution of microstates. In particular, for temperatures far below the
interaction energy, $U/t$ sensitively determines the number of single occupied
sites. The value of $U/t$ reported in the main text comes from a best fit of
the experimental data, which is in reasonable agreement with an independent
spectroscopic calibration [$U/t=4.4(5)$].

To prepare the $2\times 2$ plaquette ground state, we initialize the ground
state of sites 1 and 4 with doublons, as pictured in Fig. 4 of the main text.
To achieve a half-filled system, we turn on the opposing diagonal sites in a
two-step sequence. First, the intensities of sites 2 and 3 are increased in 5
ms with an exponential time constant of 1 ms, until they are at $\sim$90%
intensity of sites 1 and 4. Then, the intensities are further increased in a
50 ms exponential ramp with a time constant of 5 ms, until the offsets are
zeroed. During this time, the magnetic field is ramped from a nearly non-
interacting system at 573 G to its final value of 631 G using a two-point
cubic spline function. A hold time of 4 ms is applied for further
equilibration. Finally, tunneling is frozen by rapidly increasing the tweezer
intensity to 2.5 times the science depth, and imaging commences.

In the main text of the paper, for the $2\times 2$ tunneling data we average
over states obtained by interchanging $|{\uparrow}\rangle$ and
$|{\downarrow}\rangle$. This is because all potentials in the system are state
independent. This reduces the statistical error due to the limited number of
experimental shots (673 shots). Since we do not equalize the disorder entirely
to zero, we do not assume any reflection symmetries that could be utilized in
the ideal system to further reduce the statistical errors. We show the full
microstate distribution without averaging in Fig. S3, and the full labelling
of microstates in Tab. 1.

Figure S3:  Full microstate data for Fig. 4 of the main text without averaging to take advantage of the spin symmetry of the Hamiltonian. Index |  $|{\uparrow}\rangle$ population |  $|{\downarrow}\rangle$ population
---|---|---
0 | 1100 | 1100
1 | 1010 | 1100
2 | 1001 | 1100
3 | 0110 | 1100
4 | 0101 | 1100
5 | 0011 | 1100
6 | 1100 | 1010
7 | 1010 | 1010
8 | 1001 | 1010
9 | 0110 | 1010
10 | 0101 | 1010
11 | 0011 | 1010
12 | 1100 | 1001
13 | 1010 | 1001
14 | 1001 | 1001
15 | 0110 | 1001
16 | 0101 | 1001
17 | 0011 | 1001
18 | 1100 | 0110
19 | 1010 | 0110
20 | 1001 | 0110
21 | 0110 | 0110
22 | 0101 | 0110
23 | 0011 | 0110
24 | 1100 | 0101
25 | 1010 | 0101
26 | 1001 | 0101
27 | 0110 | 0101
28 | 0101 | 0101
29 | 0011 | 0101
30 | 1100 | 0011
31 | 1010 | 0011
32 | 1001 | 0011
33 | 0110 | 0011
34 | 0101 | 0011
35 | 0011 | 0011
Table 1:  The basis microstates for the data reported in Fig. 4 of the main
text. The order of the sites used corresponds to the numbering in Fig. 4(a).

##  5 Entropy and Temperature

In order to fit the temperature of our system, we use the canonical ensemble
as post-selection enables us to remove number fluctuations. Thus, with free
parameters temperature $T$ and $U/t$, the partition function is the sum over
all energy eigenstates $\Omega$, with

$Z\left(U/t,T\right)=\sum_{\Omega}\exp\left(\frac{-E_{\Omega}(U/t)}{k_{b}T}\right).$
(S6)

For each temperature and interaction, we create an expected distribution of
microstates, and we compare to the experimental data using a weighted least
squares fit. We find that the weighted sum of the squared residuals reaches a
local minimum (which depends very weakly on temperature) at $U/t=3.4(2)$,
where the errorbar is extracted from 500 bootstrapped samples of the data. We
find that the temperature fit loses sensitivity below $k_{\mathrm{B}}$$T\sim
0.3\bar{t}$ (Fig S4). This corresponds to an entropy range of [0,0.09]
$k_{\mathrm{B}}$ per particle. Exact diagonalization performed by in the
Python package Quspin [49] suggests that for our ramp parameters, the entropy
gain would be 0.02 $k_{\mathrm{B}}$ per particle.

Figure S4:  Entropy and least squares residual as a function of temperature.
Below around $k_{\mathrm{B}}$$T=0.3\bar{t}$, the least squares residual is
minimized, corresponding to an entropy range of [0,0.09] $k_{\mathrm{B}}$ per
particle.

In future studies, we wish to implement larger 2D ladder systems, preparing
ground states of half-filled systems by the aforementioned technique of
adiabatically ramping on additional tweezer sites. We use exact
diagonalization on numerically tractable systems to estimate the many-body
gaps sizes and therefore the feasibility of this approach. For instance, on a
$2\times 5$ triangular lattice at half filling, at $U/t=6$, we can perform an
adiabatic ramp in 100 ms that will only lead to an increase of 0.06
$k_{\mathrm{B}}$ per particle. For a $2\times 6$ lattice with the same
parameters, the increase is 0.07 $k_{\mathrm{B}}$. These values were obtained
by calculating dynamics with experimentally realistic ramp parameters. We are
unable to perform exact diagonalization with our computing resources for
larger 2D systems.