Non-intuitive Computational Optimization of Illumination Patterns for Maximum Optical Force and Torque

TL;DR

This paper presents a numerical framework for designing structured illumination to significantly enhance optical force and torque on micro- and nano-scale objects, outperforming traditional methods.

Contribution

It introduces a novel computational approach combining Bessel basis, boundary element method, and local optimization for optimized illumination design.

Findings

Achieved 20-fold increase in optical torque compared to circular polarization.

Demonstrated effectiveness across 2000 random initial configurations.

Analyzed the impact of multipolar plasmon resonances on force and torque enhancement.

Abstract

This paper aims to maximize optical force and torque on arbitrary micro- and nano-scale objects using numerically optimized structured illumination. By developing a numerical framework for computer-automated design of 3d vector-field illumination, we demonstrate a 20-fold enhancement in optical torque per intensity over circularly polarized plane wave on a model plasmonic particle. The nonconvex optimization is efficiently performed by combining a compact cylindrical Bessel basis representation with a fast boundary element method and a standard derivative-free, local optimization algorithm. We analyze the optimization results for 2000 random initial configurations, discuss the tradeoff between robustness and enhancement, and compare the different effects of multipolar plasmon resonances on enhancing force and torque. All results are obtained using open-source computational software available online.

Figures

Equations

$${\psi }_m(r,\phi ,z)=J_m(k_{t}r)\exp (im\phi +i{k}_{z}z),$$

$${\mathbf{M}}_i$$

$${\mathbf{N}}_i$$

$${\mathbf{E}}_{\text{inc}}(r,\varphi ,z)=\sum _{i=0}^N{a}_i{\mathbf{M}}_i+b_i{\mathbf{N}}_i,$$

$$\underset{\begin{array}{c}\text{fixed}\\ \text{BEM matrix}\end{array}}{\underbrace{\mathbf{M}}}\underset{\begin{array}{c}\text{output}\\ \text{current}\end{array}}{\underbrace{\mathbf{c}}}=\underset{\begin{array}{c}\text{rapidly updated}\\ \text{input field}\end{array}}{\underbrace{\mathbf{f}}}$$

$${\text{FOM}}_F$$

$${\text{FOM}}_T$$

$${\mathbf{E}}_w=\sum _{i=0}^m(a_i+{\delta }_{ai}){\mathbf{M}}_i+(b_i+{\delta }_{bi}){\mathbf{N}}_i,\left(|\delta |\le w\cdot |a,b{|}_{\infty }\right),$$

Full text

1Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

2Department of Mathematics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

2Currently with the Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USA

*nicfang@mit.edu

OSCIS code (350.4855) Optical Tweezers or optical manipulation;(140.7010); (090.1760) Computer holography; Diffraction and Gratings: Optical Vortices; (250.5403) Plasmonics; (290.2200) Extinction

1 Introduction

We show how large-scale computational optimization [1, 2, 3] can be used to design superior and non-intuitive structured illumination patterns that achieve 20-fold enhancements (for fixed incident-field intensity) of the optical torque on sub-micron particles, demonstrating the utility of an optimal design approach for the many nanoscience applications that rely on optical actuation of nanoparticles [4, 5, 6].

Recent advances in nanoparticle engineering [7, 8] and holographic beam-generation via spatial light modulators (SLMs) [9, 10, 11, 12] and other phase-manipulation techniques[13, 14, 15, 16] have created many new degrees of freedom for engineering light–particle interactions beyond traditional optical tweezers. Enhanced and unusual optical forces and torques can be engineered by designing material objects [17, 18, 19, 20, 21] and/or structured illumination, with the latter including “tractor beams” [22, 23] and beams carrying optical angular momentum [24, 25, 26, 27]. These increased degrees of freedom pose an interesting design challenge: for a given target object, what is the optimal illumination pattern to produce the strongest optical force or torque? While a small number of Gaussian beam parameters can be manually calibrated for optimal performance using manual trial-and-error [28], an arbitrary 3D vector field requires a more targeted approach. Moreover, optimization of a 3d vector field is highly nonconvex by nature, and posesses many local optima due to wave interference and resonance. When exploring so many parameters, a large number of scattering problems must be solved efficiently, which requires careful design of the optimization framework.

Research in computational optimization of optical actuation has focused on the design of new material geometries [29, 30] and on the improvement of multiplexed optical traps for microscale dielectric particles (holographic optical tweezers) [31, 32, 33, 34, 35, 36, 37]. However, no computational method has been available to design structured illumination for unconventional target objects that are nonspherical, lossy, or nanometer-scale, thereby requiring a costly full-wave numerical simulation for computing optical force and torque.

We present a compact and rapid computational framework to optimize structured illumination for the mechanical actuation of an arbitrary target object. We combine (i) a compact Bessel-basis representation (Sec.2.1); (ii) a numerical solver based on boundary element method (BEM) that discretizes the surfaces of a 3d scattering problem to form a BEM matrix, and solves hundreds of thousands of incident-field configurations using the same matrix (Sec.2.2); (iii) an appropriate optimization algorithm that exploits the smoothness of the nonconvex and nonlinear optimization problem (Sec.2.3); and (iv) a suitable figure of merit (FOM) and optimization constraints (Sec.2.4). As a result, we rapidly attain many-fold improvements in optical forces and torques over random field or plane wave illuminations (Sec.3.1), and discuss the tradeoff between enhancement and robustness of optimization (Sec.3.3). Furthermore, a given material object may have scattering resonances at various frequencies, and the choice of frequency for the incident field has several important implications. When comparing interactions with two different resonances, we were able to distinguish the impact of the change in the resonant field pattern from the change in the resonance lifetime. Controlling for the change in lifetime, we found that torque seems to favor higher-order (e.g. quadrupole) resonances with greater angular momentum, while force seems to favor lower-order (e.g. dipole) resonance with greater field intensity within the particle (Sec.3.2).

2 Optimization Framework

A structured-illumination optimization aims to find the best 3d vector field, i.e., the one that maximizes a desired FOM for a given scattering problem (see Fig.1). Our choice of the Bessel basis expansion is described in Sec.2.1. Sec.2.2 and Sec.2.3 discuss the BEM solver and the optimization algorithm. Lastly, Sec. 2.4 explains force and torque FOMs and the corresponding optimization constraints. The entire numerical framework is implemented with C++.

2.1 Analytical Representation of Structured Illumination

The computational design of structured illumination requires a compact analytical representation of an arbitrary 3d vector field. The vector field will contain spatial variations in both intensity and phase, and must satisfy the vector wave equation [38]. The electric field $\mathbf{E}$ can be represented using a basis expansion $\mathbf{E}=\sum_{i=0}^{N}c_{i}\boldsymbol{\phi}_{i}$, where the complex scalar coefficient $c_{i}$ determines the relative intensity and phase of each mode $\boldsymbol{\phi}_{i}$.

The choice of coordinates and the basis functions $\boldsymbol{\phi}_{i}$ depends on the problem geometry. While the spherical coordinate system is a common choice in Mie scattering [39, 40], it requires a very large number of modes $N$ to represent light propagating along a linear axis and potentially interacting with flat substrates or SLMs. The cartesian coordinate system is also ill-suited because it requires large $N$ to describe laser beams with a finite radius. We find that the cylindrical coordinate system [41, 42] is well suited, requiring a small number of modes to describe structured illumination with varying distributions of linear and angular momentum.

Among the wide menu of cylindrical basis functions (e.g., Bessel, Laguerre-Gaussian, Hermite-Gaussian, and so forth) [42], we choose the Bessel basis [43, 44] for its compact analytical expression, derived from the scalar generating function

[TABLE]

where $J_{m}$ is the mth-order Bessel function, $k_{t}$ is the transverse wavevector in $\hat{r}$, and $k_{z}$ is the longitudinal wavevector in $\hat{z}$, satisfying $k_{t}^{2}+k_{z}^{2}=2\pi/\lambda$. The ratio between $k_{t}$ and $k_{z}$ specificies the numerical aperture of the basis, $\textit{NA}=\tan^{-1}(k_{t}/k_{z})$. Higher NA represents greater transverse momentum that can increase optical torque. But the permissible range of NA is often dictated by experimental considerations, and the range of NA in the optimization can be set accordingly.

Taking spatial derivaties of Eq.(1) gives

[TABLE]

where $\mathbf{u}_{z}$ is the unit vector in $\hat{z}$, and $\mathbf{\mathbf{M}}_{i}$ and $\mathbf{N}_{i}$ are the $i$th bases for azimuthal and radial polarizations, respectively (right inset of Fig.1). The incident electric field $\mathbf{E}_{\text{inc}}$ can be expressed as:

[TABLE]

where $a_{i}$ and $b_{i}$ are the complex scalar coefficients. The Bessel basis produces the most compact expressions for $\mathbf{\mathbf{M}}_{i}$ and $\mathbf{N}_{i}$ because the magnitude of $\psi$ does not vary with $z$, reducing $\partial/\partial z$ terms in Eqs.(2,3).

Note that we intentionally decouple our optimization framework from the idiosyncratic differences in the spatial resolution of the SLMs. A wide variety of experimental methods (e.g., superposed pitch-fork holograms) [9] can be used to generate beams expressed as Eq.(4), for a finite $N$ and NA. In this paper, we consider numerical apertures with opening angles $\leq 10^{\circ}$ and 12 basis functions ($N=5$).

2.2 Numerical Solver

The optimization process itself has no restrictions on the choice of the numerical solver, so the biggest consideration is the computational cost: the smaller the the better. We choose the Boundary Element Method (BEM) [45, 46] for several reasons. In comparison to other scattering methodologies such as the finite-difference or finite-element methods, BEM is particularly well-suited to the type of large-scale optimization problem requiring a rapid update of the incident field for a given geometry and a given wavelength $\lambda_{\text{opt}}$.

[TABLE]

In Eq.(5), the BEM matrix $\boldsymbol{M}$ remains fixed for a given geometry and frequency, while the column $\mathbf{f}$ representing the incident field is rapidly updated at each step of the optimization process. This allows hundreds of thousands of scattering configurations to be computed on the order of a few hours. In addition, BEM projects the 3d scattering problem onto a 2d surface mesh, thereby reducing the computation volume by a factor of 1/2000 in the nanoparticle scattering problem we consider. Lastly, recent improvements [47] have significantly increased the speed with which optical froce and torque can be computed in BEM.

2.3 Optimization Algorithm

Structured-illumination optimization is nonlinear and nonconvex, such that searching for a global optimum is prohibitively expensive. Therefore we choose a local algorithm with random starting points. We choose one of the simplest solutions available: constrained optimization by linear approximation (COBLYA) [2], a derivative-free algorithm that exploits the smoothness of the problem. An open-source implementation of COBYLA is available through NLopt [3].

2.4 Figure of Merit and Optimization Constraints for Optical Force and

Torque

We consider two types of optical actuation with respect to the object coordinate (left inset of Fig.1); the force $F_{z}$ and torque $T_{z}$. In order to avoid the optimizer from increasing the brightness of the beam indefinitely, we choose to divide the force and torque by the average incident-field intensity on the particle surface ($I_{\text{avg}}=\overline{|E_{\text{inc}}|^{2}}/2Z_{0}$, where $Z_{0}$ is the impedance of free space), which is easily computed in BEM. We choose the incident-field rather than the total-field intensity to avoid penalizing high extinction efficiency. $I_{\text{avg}}$ is measured on the particle surface because we want to account for the portion of the beam that interacts with the target particle, rather than the entire beam. We choose nondimensionalized figures of merit:

[TABLE]

where the constants in parentheses reflect ideal single-channel scattering. The largest [48, 49, 50] scattering cross-section into a single (spherical harmonic) channel is $3\lambda^{2}/2\pi$, which when multiplied by single-photon changes in linear ($2\hbar k$) and angular ($\hbar$) momentum per photon, divided by the photon energy ($\hbar\omega$), yields the constants in Eqs.(6,7).

Optimization constraints can be added to suppress actuation in undesired directions. We suppress actuations in directions other than $\hat{z}$ using smooth constraints: $\left(|\mathbf{F}|^{2}-F_{z}^{2}\right)/|\mathbf{F}|^{2}\leq 0.01$ and $\left(|\mathbf{T}|^{2}-T_{z}^{2}\right)/|\mathbf{T}|^{2}\leq 0.01$, where the limiting value 0.01 is set to ensure that $F_{z}$ and $T_{z}$ exceed $99\%$ of the force and torque magnitudes $|\mathbf{F}|$ and $|\mathbf{T}|$, respectively.

3 Results and Discussion

We demonstrate our illumination-field optimization framework on the gold nanotriangle illustrated in Fig.1. Our previous work [51] analyzes the optical force and torque on such a particle for circularly-polarized (CP) planewave illumination. CP planewave is a common incident-field choice [52, 53, 19, 20] for torque generation due to its intrinsic spin angular momentum, but we find in our computational optimization that highly optimized field patterns can show 20x improvement of $\text{FOM}_{T}$. The wavelength of illumination in each optimization, $\lambda_{\text{opt}}$, is chosen to correspond to the plasmonic resonance wavelengths of the model particle.

Sec.3.1 presents the distribution of 2000 local-optimization results and discusses the optimized field-patterns. Sec.3.2 analyzes the wavelength-dependence of optical force and torque for optimized illuminations, based on the choice of $\lambda_{\text{opt}}$, and compares the results with the reference force and torque from CP planewave. Lastly, the tradeoff between robustness and enhancement is discussed in Sec.3.3.

3.1 Illumination-field Optimization from 2000 Randomly Selected Initial

Configurations

The illumination-field design space is nonconvex and littered with local optima, due primarily to wave-optical interference effects. We survey this broader design space by restarting our local-optimization algorithm 2000 times with randomly selected initial configurations that are constructed using Eq.(4), where the complex coefficients $a_{i}$ and $b_{i}$ are uniform random numbers bounded by $|a_{i}|,|b_{i}|\leq 1$. The results are summarized in Figs.2-3 at $\lambda_{\text{opt}}=1028\text{nm}$ (dipole resonance) and 625nm (quadrupole resonance), respectively. At both wavelengths, we find that more than 50% of local optimizations from random starting points can achieve over 5x enhancement of $\text{FOM}_{T}$ compared to CP planewave reference, and that the optimized field patterns contain various combinations of Bessel-basis modes without a systematic convergence to one over the others. The distributions are plotted in log-scale to increase the visibility of small bins.

In Fig.2, the median $\text{FOM}_{T}$ at 0.913 is very close to the best $\text{FOM}_{T}$ at 1.01 and the distribution is predominantly concentrated to the right: the 4 rightmost bars represent 69% of all samples, which all acheive over 5x enhancement of $\text{FOM}_{T}$ compared to CP planewave reference at 0.169 (marked with a red triangle). The insets show the optimization results from two different starting points – random field (top) and CP planewave (bottom) – that reach similar optimized field patterns. We also observe that a variety of other patterns, dominated by different combinations of Bessel-basis modes, can produce a nearly identical or superior FOM.

In Fig.3, the final $\text{FOM}_{T}$ distribution is more dispersed between the median at 3.934 and the best at 10.94, which respectively achieve over 5x and 14x enhancement compared to CP planewave reference at 0.779. As in Fig.2, the results concentrate heavily around the median; however, in Fig.3 a small number of samples achieve a remarkable improvement above 14-fold. The top inset shows four different field patterns that produce a nearly identical $\text{FOM}_{T}$ above the median, and the bottom inset shows the field-pattern with the highest $\text{FOM}_{T}$. A comparison of all optimized field patterns at 1028nm and 625nm shows that the latter contains more higher-order Bessel-basis contributions.

3.2 Dependence on Illumination Wavelength

We further investigate the influence of $\lambda_{\text{opt}}$ by plotting optical force and torque per incident-field intensity $I_{\text{avg}}$ as a function of illumination wavelength. In Figs.4A-4B, the reference planewave force spectrum (black dashed line) is identical in both plots, clearly dominated by a broad dipole resonance with smaller peaks at higher-order resonances. Through illumination-field optimization, we can enhance the force at the dipole mode while suppressing higher-order modes (Fig.4A) and also enhance the force at the quadrupole mode while suppressing the dipole mode (Fig.4B).

In Figs.5A-5C, the reference planewave torque spectrum (black dashed line) is identical in all three plots and peaks at both dipole and quadrupole resonances with nearly equal heights (explained in detail in [51]). Illumination-field optimization at dipole resonance and off-resonance achieve a similar 6x-boost at the chosen $\lambda_{\text{opt}}$ value without suppressing the quadrupole resonance. The optimized total fields in Fig.5D at 1028nm and 805nm both exhibit a $4\pi$ phase change around the circumference of the particle, resembling a quadrupole resonance.

In 5C, on the other hand, the best optimization at quadrupole resonance achieves a remarkable 20x improvement while suppressing much of the dipole resonance; the median optimization also achieves 12x improvement while suppressing the dipole resonance to a lesser extent. The optimized total field in Fig.5D at 625nm shows a distinct, highly resonant distribution.

When comparing interactions with two different resonances, we were able to distinguish the impact of the change in resonant field pattern from the change in the resonance lifetime. Controlling for the change in lifetime, we found that torque seems to favor higher-order (e.g. quadrupole) resonances with greater angular momentum, while force ($F_{z}$), which closely correlates with extinction power, seems to favor lower-order (e.g. dipole) resonance with greater field intensity within the particle. With the use of structured illumination, higher-order resonances can be excited more effectively, which contributes to higher optical torque after optimization.

Our numerical optimization framework allows a systematic search of the illumination-field design space to maximize force and torque on lossy, non-spherical particles with multipolar scattering channels. In addition, we think a rigorous analytical study of the fundamental upper bounds on opto-mechanical responses, similar to the analysis performed on light extinction [54, 55], would be useful in the future.

3.3 Robustness of Optimization

Experimental generation of the designed illumination via SLMs may suffer various types of manufacturing errors [56]. Fig.6 shows the tradeoff between enhancement and robustness to experimental errors. The fractional error is added to the beam using

[TABLE]

where $\delta$ is a complex-valued random error bounded by the magnitude of the largest coefficient multiplied by the fractional weight $0\leq w\leq 1$. Fig.6 shows the decrease in $\text{FOM}_{T}$ as a function of $w$. At $\lambda_{\text{dip}}$, increasing $w$ from $1\%$ to $10\%$ decreased the best FOM from $0.994$ to $0.517$, and the median FOM followed a similar trend. On the other hand, $\text{FOM}_{T}$ of the best optimized field at $\lambda_{\text{quad}}$ dropped from $8.624$ to $1.433$, and the median optimized field changed from $3.92$ to $3.237$. Fig.6 demonstrates a clear tradeoff between field enhancement and error tolerance, as one might expect due to the need to couple strongly to the underlying particle resonances. For the tolerance requirements of a given experimental setup, the approach we outline here could easily be adapted to a robust-optimization framework [57, 58] in which the expected variability is included and optimized against.

4 Conclusion

We present a numerical framework for computer optimization of structured illumination that maximizes optical force and torque on arbitrary scatterers, and show a 20-fold enhancement in optical torque per intensity on an example plasmonic nanoparticle, compared to a circularly polarized planewave. Previously, the major bottleneck has been the cumbersome computation. We overcome this bottleneck with a compact cylindrical Bessel basis and a fast boundary element method. We are optimistic that such computational framework for 3d vector fields can be generalized and applied to other design problems in opto-mechanics, nanophotonics, and 3d imaging.

Acknowledgement

The authors thank George Barbastathis for helpful discussions.