Ultrafast and high resolution spatial light modulation for cold atoms

In this work we define the frame rate as the inverse of the rise time. We define the spatial resolution of the system in terms of the $\frac{1}{e^{2}}$ radii given by an elliptical Gaussian fit of the PSF shown in Fig. 1(d).

The optical design of the system (Fig. 5) is compatible with simultaneous operation over $>10\penalty\$THz of bandwidth. Here, as a proof of principle, our signal path is composed of a single Mach-Zehnder modulator with 40 GHz of bandwidth, which is paired with a $780\penalty\$nm laser source that is tunable over $>300\penalty\$GHz.

After passing through an MZM, the beam is expanded along one axis in a 4:1 anamorphic prism pair. This ensures that the beam does not diffract significantly in the $x$ axis (the axis orthogonal to the VIPA deflection) as it traverses the VIPA. Note that diffraction in the $y$ axis (along which the VIPA deflects) does not substantially modify the resolution of the VIPA, as discussed in the main text.

The VIPA has a gradient coating (Fig. 6), resulting in an approximately uniform intensity envelope across all emitters. The output of the VIPA is expanded in the $y$ axis using a 4:1 cylindrical telescope. This magnification is chosen to tune the effective resolution of the grating, and produces an approximately square beam in the 29th order of the grating.

The resulting output has significant aberrations, primarily due to the many passes through the VIPA. By far the largest aberration is associated with a slight wedge of the VIPA, which results in a departure angle for each emitter that scales linearly with the number of passes through the VIPA. In the direction of the VIPA deflection, this results in a chirp in the positions of the emitters, and a slight nonlinearity that can be corrected with the appropriate conversion between frequency and position in the signal path. In the orthogonal direction, this results in astigmatism, which can be corrected by rotating one of the cylindrical lenses in the system.

The output of the grating is imaged onto an LCOS SLM via a 4f telescope with a demagnification of $1/4$, which allows one to compensate for the remaining aberrations in the system. The SLM can additionally be used for amplitude modulation, to produce multiple copies of the image for parallelized operation across multiple logical qubits (e.g. for magic state distillation), and to scan the image produced by the system across a fiber for the fast video recording scheme described in Appendix J.

In atom trapping applications, the output of the SLM could be imaged onto the atoms using a high NA microscope objective. Here, the output is imaged onto a camera for characterization.

The overall efficiency of the system is 5.8%, which is primarily limited by geometric constraints and the selected grating. An optimized design with a custom VIPA and grating could achieve efficiencies of $>50\%$.

To demonstrate the capabilities of our approach unconstrained by on-hand hardware, we apply several post-processing steps. These are all equivalent to operations that can readily be performed optically in a future setup. First, in Fig. 1(c) we measure the optical transfer function by scanning the laser frequency and recording the PSF across the field of view. We then stitch together the PSFs with adjusted weighting to produce an arbitrary image. Similarly, in Fig. 3(a) we modulate individual FSRs and then compute the resulting image by summing the individual photos together. These are both equivalent to modulating a wider RF bandwidth or wavelength division multiplexing. Second, when modulating, we crop the carrier and negative sidebands from photos. This is equivalent to optically filtering or single sideband modulation with an extinction ratio surpassing our available MZM. Finally, we compensate for ellipticity of the PSF by stretching the $y$ axis by $w_{x}/w_{y}$, based on a PSF measurement near the center of the field of view under the same conditions. This is identical to a cylindrical telescope. For each figure we fit a representative PSF under current conditions to determine the appropriate scaling factors. An example fit can be seen in Fig. 8(c). These parameters are listed in Table 1 where $w_{x}$ and $w_{y}$ are the $\frac{1}{e^{2}}$ radii closest to the $x$ and $y$ axes respectively, and $\theta_{\text{psf}}$ is the tilt of the PSF principal axes relative to the camera axes.

Note that these elliptical Gaussian fits, as well as those in Fig. 10, are performed by first fitting or estimating the PSF, then refitting restricted to a region two waists in radius around the predicted spot to focus on the central lobe. This slightly reduces the measured waist.

In this work we use two different signal paths optimized for different goals. Setup (A) uses a Tektronix AWG70001A AWG operating at 48.8 GS/s with a $20\text{\,}\mathrm{GHz}$ bandwidth and 7-bit vertical resolution. This is paired with two Lotus Systems LNA2G18G low noise amplifiers (LNAs). This setup operates between $2$ GHz and $18$ GHz with a theoretical maximum output power of $24$ dBm, and is used for characterizing large bandwidth operations such as moving spots and measuring modulation bandwidth, but struggles with linearity and vertical resolution. Vertical resolution is critical in our system for intensity control, since the finite resolution is shared across all tweezers. Likewise, crosstalk due to nonlinearity limits the effectiveness of our calibration procedure which assumes independence between tones.

Setup (B) uses a Xilinx ZCU216 RFSoC development board running the QICK firmware Ding et al. 2024. This setup operates at 9.8 GS/s with 14-bit vertical resolution and is paired with a Mini-Circuits ZVE-3W-83+ $2$–$8$ GHz high-power linear amplifier and Mini-Circuits VHF-4600+ and VLF-8400+ filters to select the second Nyquist zone. This setup operates between $5.4$–$8$ GHz and has linearity limited only by the MZM itself.

To calibrate the mapping between RF power and optical intensity, we play a sum of tones corresponding to all targeted tweezers and measure the resulting intensity by fitting an elliptical Gaussian to the spot. It is helpful to crop the fit to an elliptical region around the predicted spot location, as in Appendix C, to minimize crosstalk. After measuring the mapping at several RF powers we fit a quadratic to the resulting curve, which is used for all subsequent amplitude modulation [Fig. 7(a)]. Since tweezers can be measured in parallel, it is possible to perform this calibration in $O(1)$ time for an arbitrarily large array. Calibrating once and computing the mean normalized RMS as a function of time reveals that without any particular effort to stabilize the system, the performance is maintained for at least $\sim 20$ minutes [Fig. 7(b)].

A remaining source of crosstalk in our current apparatus comes from the finite extinction ratio of the MZM and the presence of negative sidebands. The $\text{sinc}^{2}$ profile expected of the VIPA without apodization has power law tails which can be significant for tweezers that are closely spaced in the $y$ axis. For demonstration we ameliorate this issue by staggering the spots in the $y$ axis. However, this significantly limits the number of spots that can be balanced. Similar crosstalk also arises from the carrier, which is worsened by its increased intensity. We expect these issues to be preventing us from reaching error below $10^{-4}$, which would otherwise be within the vertical resolution of the RFSoC. Single sideband modulation as described in Appendix C and apodization of the PSF would resolve this issue.

To estimate the optimal SLM correction, we apply a hill climbing procedure to each coefficient of a rectangular Legendre polynomial, up to total degree $5$, trying to maximize the ratio of the intensity within a $1$-pixel-radius region to that within a $20$-pixel-radius region. After each hill climbing iteration we fit a quadratic at the optimum which we accept if it improves the score, and repeat over these terms until convergence. This procedure converges quickly and is effective at recovering diffraction-limited performance. We include a representative optimization curve and resulting phase mask in Fig. 8, along with line cuts of the resulting PSF. Future work may explore more sophisticated optimization techniques, including machine learning approaches or segment-based interferometric wavefront sensing Zuo et al. 2022; Deutsch et al. 2008.

For generating waveforms and computing spatial resolutions we use a basic model of the system’s dispersion, where we assume linear dispersion in both axes and orthogonality between the VIPA and the grating. To determine model parameters, we play a sequence of known frequencies from a DS Instruments SG22000PRO signal generator and fit to tweezer locations on the camera [Fig. 9(a)]. By measuring the dispersion when moving within a single FSR and projecting the displacement onto the vector between two adjacent orders, we can estimate the position shift due to the VIPA alone. Based on the frequency needed to move between these orders we extract the FSR, and assuming the grating is orthogonal to the VIPA, we fit its dispersion across multiple FSRs. This procedure predicts tweezer locations within 1 camera pixel [Fig. 9(b)], which is sufficient for our purposes, although it neglects nonlinearities in the VIPA and inexact orthogonality between the VIPA and the grating, so we expect systematic errors. Future megapixel-scale VIPA SLMs will likely demand more detailed, nonlinear models of the VIPA. Such models have been explored in spectroscopy literature Xiao et al. 2004 and can be directly applied to our system.

We extract relevant parameters for resolution which we tabulate in Table 2. We do not report uncertainties as shot noise in the camera is negligible and deviation is dominated by systematic errors in the model. These values together with the measured waists from Figs. 1(c) and 1(d) in Table 1 are used to estimate the number of resolvable spots in the system. In particular, the frequency step $f_{\alpha}^{\prime}$ which displaces by a waist is given by $\frac{\delta\vec{r}_{\alpha}}{\delta f_{\alpha}}f_{\alpha}^{\prime}=\vec{w}_{\alpha}^{\parallel}$ where $\vec{r}_{\alpha}$ is the displacement vector for the grating ($\alpha$ = g) or VIPA ($\alpha$ = v) and $\vec{w}_{\alpha}^{\parallel}$ is the projection of the waist along the displacement vector (computed using $\theta_{\text{psf}}$ in Table 1). Finally, the number of resolvable spots is given by $N_{\text{v}}=\left\lfloor\frac{\Delta f_{\text{v}}}{f_{\text{v}}^{\prime}}\right\rfloor$ and $N_{\text{g}}=\left\lfloor\frac{$10\text{\,}\mathrm{THz}$}{f_{\text{g}}^{\prime}}\right\rfloor,\left\lfloor\frac{$286.9\text{\,}\mathrm{GHz}$}{f_{\text{g}}^{\prime}}\right\rfloor,\left\lfloor\frac{$40\text{\,}\mathrm{GHz}$}{f_{\text{g}}^{\prime}}\right\rfloor$ for the full C+L band, field-of-view limited span, and the 40 GHz span respectively.

Fitting elliptical Gaussians to the spots shown in Fig. 1(c), we find minimal variation in the waist across the field of view [Figs. 10(a)–10(e)]. Note that these waists are measured with a static SLM correction near the center of the field of view. Some broadening is observed near the edges due to clipping of the beam on the 2-inch optics used in this work.

To create a pure translation of a tweezer in the horizontal direction, we compensate for any tilt in the grating which would otherwise cause a change in vertical position when naively moving from $f_{0}$ to $f_{0}+\Delta f_{\mathrm{vipa}}$ by adding,

Next, when doing handoffs across FSRs there is always broadening because the resulting spot is an interpolation of spatially separated tweezers. In this work, up to corrections for amplitude rolloff of the RF path with frequency, we use the functional form $V_{1}^{2}=1-s$ and $V_{2}^{2}=s$ with $s\in[0,1]$ for the two tones to linearly move the center of mass. A downside to this approach is that it results in a non-constant waist during the move. In future work, with three overlapping tweezers, one can maintain constant broadening throughout (defined based on the second central moment) while linearly moving the center of mass. To achieve this with a triplet of adjacent spots $(V_{k-1},V_{k},V_{k+1})$ we can parametrize based on the offset of the center of mass $\Delta=x_{\mathrm{COM}}-x_{k}\in[-d/2,d/2]$ and the distance between adjacent spots $d$,

This gives linear motion of the COM and a constant waist of $w=\sqrt{w_{0}^{2}+d^{2}}$ during the move. At the start/end of the move a two spot handoff can be used. Note that this is why there is no fundamental broadening when chirping within an FSR as $d\approx 0$.

This result implies a tradeoff between the bandwidth needed to cover a given distance along the grating axis and the resulting waist size during the move. However, this scaling is favorable in the sense that we must cover $N_{w}=\frac{w_{0}}{d}$ FSRs to move by one waist, so for small $\frac{d}{w_{0}}$ we get broadening,

So broadening reduces quadratically with bandwidth devoted to the handoff.

To record the dynamics of the dSLM we use the LCOS SLM already present in our system to scan the image plane across a $0.22$ NA bare fiber, whose output is focused onto a DC-400 MHz Thorlabs APD430A. This is referred to as a single-pixel camera Gibson et al. 2020. We average each trace over $1000$ experiments per pixel to determine a mean and uncertainty for the intensity at each timestep.

With the optically conjugate CCD camera we reference the fiber location to a pixel by maximizing fiber coupling of a single spot and recording its location on the camera, doing a Gaussian fit for subpixel accuracy. With this information, we can compensate for drifts in beam pointing and verify the point we are imaging is aligned within $0.2$px to the fiber location at each pixel we record.

Traces are captured on a $2\text{\,}\mathrm{GHz}$ LeCroy WaveRunner 204 MXi-A with a $25$ sample linear phase finite impulse response (FIR) filter of $290$ MHz bandwidth. In total the detector chain has a net bandwidth of $(($290\text{\,}\mathrm{MHz}$)^{-2}+($400\text{\,}\mathrm{MHz}$)^{-2}+($2\text{\,}\mathrm{GHz}$)^{-2})^{-1/2}\simeq$233\text{\,}\mathrm{MHz}$$. This implies a detector rise time of $1.5\text{\,}\mathrm{ns}$, meaning we overestimate the intrinsic dSLM rise time by only $\simeq 0.1$ ns.

Although the frame rate is directly defined by the rise time, a related question is what spectral content can be transmitted through the system, for example to modulate links for the Hubbard scheme in Fig. 3. In the dSLM, we can consider a sinusoidal amplitude modulation in the steady state as the sum of a carrier and two sidebands. As the frequency increases, the two sidebands are increasingly dispersed, effectively attenuating the modulation.

To quantify this, we apply a sinusoidal modulation with a normalized amplitude from $0.2$ to $1.0$ on top of a carrier frequency of $9.5$ GHz for the maximum number of periods that can fit in $1.01$ $\mu$s. For each modulation frequency, this is followed by a pulse to amplitude $1.0$ for normalization, which is also used for the rise time measurement. We record near the center of the resulting spot on our single-pixel camera and accumulate traces at frequencies between $1$ MHz and $50$ MHz, and observe the rolloff in amplitude with increasing frequency (Fig. 11).

Although the response time of the system is set by the time for light to traverse the VIPA, dynamics beyond this limit can be observed by considering the transient response of light propagating through the system. An example of this is shown in Fig. 12, where we apply a $\pi$ phase jump to a tweezer and observe the time evolution in the image plane. At the time when half of the light with the phase jump has traversed the VIPA, the resulting pattern will be the interference of two spots with halved resolution, producing a dark region at the center of the original tweezer. This is desirable when one wants to momentarily extinguish a tweezer faster than would be possible with the regular bandwidth, for example to avoid light shifts during a gate operation. More sophisticated control of the phase and amplitude of the light in the VIPA can allow for more complex dynamic operations in future work.

Beam propagation through the VIPA is simulated using angular-spectrum propagation. To identify the appropriate incoupling beam parameters (namely beam waist, position, and angle) a disorder-free 1D simulation is used, allowing for optimization via gradient descent or brute force search. The high reflection cutoff on the input facet is modeled as a hyperbolic tangent function in reflectivity with a transition edge of $2\penalty\ \mu m$, as is representative of high quality commercial VIPAs. To certify performance in the presence of disorder, a 2D simulation is used, where surface roughness is modeled as a spatially correlated displacement map on each surface. To model a static correction of the resulting phase errors with an LCOS SLM, we sample the phase errors at a single reference frequency, and apply a correction that is discretized to a $4160\times 2464$ pixel grid (representative of the highest resolution LCOS SLMs that are readily available). Although the SLM pattern remains fixed, as the laser frequency is varied both the phase error from the VIPA and the correction applied by the SLM shift slightly. For very large frequency changes, the correction applied by the SLM will no longer be appropriate.

To simulate tweezer arrays generated by the dSLM we assume that each tweezer has a Gaussian beam profile of equal width $w_{0}$. Position $\mathbf{r}_{i}$ and depth $D_{i}$ are free parameters that are controlled by the RF drive. Because neighboring tweezers are well separated in frequency in the dSLM setup, we assume an incoherent sum of intensities of the tweezer traps. The total potential is therefore:

We adapted the discrete variable representation (DVR) method used in Refs. Wall et al. 2015; Wei et al. 2024 to construct the eigenbasis of the tweezer array. We assume a 2D tweezer potential with fixed confinement in the $z$ direction across all sites. The DVR discretizes the 2D plane and solves the Schrödinger equation given the optical potential. Figs. 13(a) and 13(b) show a convergence test in a square lattice geometry of 15-by-15 tweezers. Unless otherwise stated, we use a square system with edge length $L=20w_{0}$, sampled on a grid with $N_{x}=N_{y}=150$ points.

For most of the geometries, we truncate the number of eigenstates to the number of tweezer traps to focus on only the lowest band. Figs. 13(c) and 13(d) show an example of the potential and energy spectrum of eigenstates in the DVR basis. The lowest band, spanned by the number of tweezers, is well-separated from the higher bands. Using the single-particle eigenstates of the optical potential for a finite tweezer array, we construct a set of maximally localized single-particle orbitals (Wannier-like functions) associated with individual tweezers. Starting from an orthonormal basis $\{\phi_{k}(\mathbf{r})\}_{k=1}^{M}$ spanning the relevant low-energy subspace, we generate a new orthonormal set $\{W_{n}(\mathbf{r})\}$ via a unitary rotation,

We determine U using the Foster–Boys localization criterion, which chooses the rotation to minimize the spatial extent of the orbitals. Defining the center of each orbital as

we maximize $F$, which is equivalent to minimizing the total quadratic spread:

Finally, we assign each localized orbital $W_{n}(\mathbf{r})$ to a specific tweezer site by matching its center $\mathbf{r}_{n}$ to the nearest tweezer position. We calculate the Hubbard parameters, including tunneling strength $t_{ij}$, on-site interaction $U_{0}$, and chemical potential $\mu$, by integrating over the real space coordinates. Figs. 13(e) and 13(f) show an example of the computed maximally localized single-particle orbitals and Hubbard parameters.

Using the DVR method, we construct a basis of maximally localized single-particle orbitals and evaluate the corresponding Hubbard parameters. Throughout this work, we use ${}^{87}\mathrm{Rb}$ as a representative atomic species. We assume a diffraction-limited optical system with center wavelength $\lambda=1064\,\mathrm{nm}$ and numerical aperture of $\mathrm{NA}=0.55$, giving a diffraction-limited beam waist:

The simulated dSLM-projected potentials are restricted to the transverse $x$-$y$ plane. Confinement in the out-of-plane direction is assumed to be harmonic, with a trap frequency of $\omega_{z}=2\pi\times 30\,\mathrm{kHz}$. The finite spatial extent of the simulated lattice leads to appreciable boundary effects: edge sites generally acquire large local energy offsets relative to bulk sites. To avoid contamination of the extracted bulk parameters, we exclude Hubbard parameters associated with the outermost sites from all subsequent analysis. These boundary sites are retained in the numerical calculation as “ghost” sites, which suppress finite-size artifacts and improve convergence of the bulk Wannier-like orbitals and Hubbard matrix elements Wei et al. 2024. Figure 14(a) compares the DVR-extracted nearest-neighbor tunneling matrix element and on-site interaction energy as functions of the base tweezer depth with the corresponding analytical expressions for a simple interfering optical lattice Bloch et al. 2008.

For the programmable square lattice, we supplement the base tweezer array with two classes of auxiliary control tweezers: $t$-control tweezers and $U$-control tweezers. These auxiliary tweezers introduce controlled perturbations to the base potential and enable local modification of the effective Hubbard parameters:

While most implementations of artificial magnetic fields modulate the local chemical potential, we are able to modulate the tunneling energy on each link to achieve the same effective Hamiltonian. This scheme allows for direct and independent access to a complex phase on each link.

One can use this modulation to realize an artificial magnetic field in a Landau gauge, which we show below. Note that the same derivation applies for the symmetric gauge.
The time-dependent Hamiltonian in the tight-binding model is given by:

We assume a sinusoidal drive of the $x$-links with frequency $\omega$, and a linear tilt in the $x$-direction that is either directly programmed using control over the local potential, or realized using a magnetic field gradient leading to a spatially varying Zeeman shift.

The driving frequency is resonant with the tilt, i.e., $\hbar\omega=\Delta_{tilt}$. In the high-frequency limit, the driving frequency is the highest energy scale in the Hamiltonian: $\hbar\omega\gg t^{(x)},t^{(y)},U$. The time-dependent Hamiltonian in the rotating frame is:

The time-independent effective Hamiltonian is obtained by performing the Magnus expansion to lowest order:

This effective Hamiltonian is the interacting Harper–Hofstadter model (HHM) in the Landau gauge when $\delta t=2t^{(y)}$. The flux $\Phi_{p,q}$ enclosed by the plaquette of (i,j), (i+1,j), (i, j+1), (i+1,j+1) is defined as:

$\Phi_{0}$ is a flux quantum. When $\Phi_{p,q}$ are identical for all plaquettes, we obtain a uniform artificial magnetic field.

In a conventional artificial magnetic field, the flux is set by the $\mathbf{k}$-vector of the light potential that generates the gauge field, i.e., $\phi_{ij}=\mathbf{k\cdot}\mathbf{r}_{ij}$. This leaves no flexibility in local degrees of freedom. In the dSLM setup, the phase on each link can be programmed independently, leading to a programmable flux on each plaquette. For example, we can insert a flux on only a single plaquette in the Landau gauge. Specifically, the HHM in the Landau gauge takes the form:

To change the flux on a single plaquette, we can change the phase difference of the top and bottom links:

However, this leads to a change of flux on the neighboring plaquette with opposite sign:

To maintain a uniform flux background, we can apply a $\delta\phi$ change to all other links along this line, extending to the boundary of the system. Note that the topology of a localized flux is precisely reflected by the fact that we needed to introduce a defect extending to the boundary of the system in order to generate it.

GitHub Copilot, Google Gemini 3.0 Pro and 3.1 Pro, and Anthropic Claude Sonnet 4.6 assisted in developing some subroutines in the control and measurement code. Copilot and Gemini were additionally used in developing portions of the angular-spectrum propagation code and scripts that loop over simulation runs. OpenAI ChatGPT 5.5 Thinking assisted with conceptualizing ideas, writing subroutines, and debugging in the Hubbard model simulation. Claude Fable 5, Opus 4.8, and ChatGPT 5.5 Thinking assisted in developing plotting code. All AI tools were used interactively, with the authors continuously validating the output.