High-fidelity holographic beam shaping with optimal transport and phase diversity

TL;DR

This paper tackles the computational difficulty of designing phase-only SLMs to produce a desired output beam and, conversely, the related task of reconstructing the input beam. By linking optimal transport to the Monge–Ampère equation, the authors derive an OT-based method that yields unwrapped, vortex-free phase estimates which serve as excellent initial guesses for iterative phase retrieval algorithms like Gerchberg–Saxton (GS) or Mixed-Region Amplitude Freedom (MRAF). They further develop phase-diversity imaging techniques, including an iterative Fourier transform (IFT) approach and one-shot/two-shot OT-based beam estimation, to recover both the amplitude and phase of the input beam from a small set of calibration images. The methods are computationally lightweight, parallelizable, and capable of high-fidelity beam shaping without GPU acceleration, with practical considerations such as memory scaling and multiscale extensions discussed. Collectively, the work provides a versatile toolkit for accurate, efficient beam shaping and robust input-beam characterization in SLM systems.

Abstract

A phase-only spatial light modulator (SLM) provides a powerful way to shape laser beams into arbitrary intensity patterns, but at the cost of a hard computational problem of determining an appropriate SLM phase. Here we show that optimal transport methods can generate approximate solutions to this problem that serve as excellent initializations for iterative phase retrieval algorithms, yielding vortex-free solutions with superior accuracy and efficiency. Additionally, we show that analogous algorithms can be used to measure the intensity and phase of the input beam incident upon the SLM via phase diversity imaging. These techniques furnish flexible and convenient solutions to the computational challenges of beam shaping with an SLM.

Figures (5)

• Figure 1: Model optical system. An input laser beam with intensity $g^2(\bm{\mathrm{x}})$ is reflected off an SLM with applied phase $\phi(\bm{\mathrm{x}})$, passes through a lens of focal length $f$ at distance $f$ from the SLM, and is then imaged on the output (camera) plane at distance $f$ from the lens, with output intensity $\tilde{G}^2(\bm{\mathrm{X}})$.
• Figure 2: Comparison of phases and output beams from various phase generation algorithms. All images are $128 \times 128$ pixels. (a) is the input beam intensity. (b) is the target output beam intensity. (c-f) are output intensities realized by the phase displayed immediately below. (g) is from GS initialized with a random phase, with RMS error $\epsilon = 13.9\%$ and efficiency $\eta = 99.13\%$. (h) is from OT; $\epsilon = 14.3\%$, $\eta = 99.96\%$. (i) is from GS initialized by OT; $\epsilon =2.58\%$, $\eta = 99.91\%$. (j) is MRAF initialized by OT; $\epsilon = 5.95\times10^{-16}$, $\eta = 85.15\%$. All iterative algorithms were run for 10,000 iterations. The MRAF hyperparameter was set by hand to 0.48. A centered $96\times 96$ pixel box was used as the MRAF signal region and the region for computing all efficiencies $\eta$.
• Figure 3: Output beam intensities resulting from various combinations of OT and IFT algorithms. The top row is a collection of target output beams from Ref. Pasienski:08. The second row is the output of our OT method with no further refinement. The third row is the output of the GS algorithm initialized with OT. The fourth row is the output of the MRAF algorithm with a tight signal region (following Ref. Pasienski:08) and initialized with OT. The fifth row is MRAF with signal region the entire field of view, again initialized with OT.
• Figure 4: Beam estimates using various phase diversity algorithms. The top row (a-d) is the beam modulus. The middle row (e-g) is the residual modulus, i.e. the difference between the modulus of the estimate and that of the ground truth. The bottom row (h-k) is the phase. (a,h) Ground truth beam. (b,e,i) One-shot beam estimate with diversity coefficient $\alpha=1.5$ ($\delta = 0.02$). (c,f,j) Two-shot beam estimate with diversity coefficients $\alpha_j=1.5$ and $\alpha_k = 0.1$ ($\delta = 0.005$). (d,g,k) IFT estimate with 15 diversity images, $\alpha=0.1,0.2,\dots,1.5$, and 1000 iterations ($\delta = 3.3\times 10^{-17}$). For visual comparison, a global phase has been chosen for each image such that the local phase in the center of the image is $0$.
• Figure 5: Beam estimation error metrics, (a) in absence of noise, (b) with additive noise, and (c) with Poissonian shot noise. Dotted and dashed lines indicate $\delta$ for the one-shot and two-shot algorithms, respectively. Solid lines show $\delta$ vs. the number of iterations of the IFT algorithm with different numbers $n$ of diversity images. Coefficients $\alpha$ in each case are as follows. $n=2: \alpha \in \{0.1, 1.5\}$. $n=3: \alpha\in \{0.1,0.8,1.5\}$. $n=4: \alpha \in \{0.1,0.6,1.0,1.5\}$. $n=6: \alpha \in \{0.1,0.4,0.7,0.9,1.2,1.5\}$. $n=15: \alpha \in \{0.1,0.2,\dots,1.5\}$. In (b), we approximate a 16-bit camera with up to two dark counts per pixel by adding to each diversity image pixel $G^2_{j,LM}$ a random value in the range $[0,2^{-15}\times \max_{PQ} G^2_{j,PQ}]$. In (c), we approximate shot noise for a 16-bit camera by letting each diversity image pixel value be a Poissonian random variable with mean $2^{16}\times G^2_{j,LM}/\max_{PQ}G^2_{j,PQ}$, where $G^2_{j,LM}$ is the corresponding noiseless pixel value. In all cases, $\delta$ is computed using all 15 noiseless diversity images.