Computational Inverse Design of Non-intuitive Illumination Patterns to Maximize Optical Force or Torque

This paper aims to maximize optical force or 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-ﬁeld 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 eﬃciently 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 conﬁgurations, discuss the tradeoﬀ between robustness and enhancement, and compare the diﬀerent eﬀects of multipolar plasmon resonances on enhancing force or torque. All results are obtained using open-source computational software available online.


Introduction
We show how large-scale computational optimization [1,2,3] can be used to design superior and nonintuitive 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 possesses 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 Figure 1: Schematic of a structured vector-field illumination, analytically represented with a vector Bessel basis. The right inset shows the distributions for electric field (color) and polarization (black arrows) of the vector Bessel basis. The illumination can be optimized to produce maximum mechanical force or torque on an example target particle. The gold nanotriangle in the left inset has edge length 400nm, thickness 40nm, and rounding diameter 30nm. The scanning electron microscope (SEM) image shows an experimental sample fabricated using electron-beam lithography. This demonstrates that such particles can be made, but all results presented in this paper are purely computational. 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 force and torque 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).

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 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++.

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 can be represented using a basis expansion E = N i=0 c i φ i , where the complex scalar coefficient c i determines the relative intensity and phase of each mode φ i .
The choice of coordinates and the basis functions φ 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 where J m is the mth-order Bessel function, k t is the transverse wavevector inr, and k z is the longitudinal wavevector inẑ, satisfying k 2 t + k 2 z = 2π/λ. The ratio between k t and k z specifies the numerical aperture of the basis, NA = arctan(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 derivatives of Eq. (1) gives where u z is the unit vector inẑ, and M i and N i are the ith bases for azimuthal and radial polarizations, respectively (right inset of Fig. 1). The incident electric field E inc can be expressed as: where a i and b i are the complex scalar coefficients. Bessel basis produces the most compact expressions for M i and N i because the magnitude of ψ does not vary with z, reducing ∂/∂z terms in Eqs. (2) and (3). We can use Eq. (4) to approximate a CP planewave configuration near the nanoparticle by setting NA= 0.01 • and a i , b i = 0, except for a 1 = 1. This approximation is used once in Fig. 2 to start the optimization from a CP planewave. All other CP planewave results are computed without approximations, using E(r, ϕ, z) = [u r + iu ϕ ] exp(iϕ + ikz) instead of Eq. (4).
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 pitchfork 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 ≤ 10 • and 12 basis functions (N = 5).

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 largescale optimization problem requiring a rapid update of the incident field for a given geometry and a given wavelength In Eq. (5), the BEM matrix M remains fixed for a given geometry and frequency, while the column 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 force and torque can be computed in BEM.

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 opensource implementation of COBYLA is available through NLopt [3].

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 avg = |E 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 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: 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λ 2 /2π, which when multiplied by single-photon changes in linear (2hk) and angular (h) momentum per photon, divided by the photon energy (hω), yields the constants in Eqs. (6) and (7). Based on the geometry and material composition of the particle, force or torque in various directions can arise even when the incident field itself does not provide an imbalance in optical momentum. Optimization constraints can be added to suppress actuation in undesired directions. We suppress actuations in directions other thanẑ using smooth constraints: |F| 2 − F 2 z /|F| 2 ≤ 0.01 and |T| 2 − T 2 z /|T| 2 ≤ 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 |F| and |T|, respectively.

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 circularlypolarized (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 FOM T . The wavelength of illumination in each optimization, λ 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 λ 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.

Illumination-field Optimization from 2000 Randomly Selected Initial Configurations
The illumination-field design space is nonconvex and littered with local optima, due primarily to waveoptical 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 | ≤ 1. The results are summarized in Figs. 2 and 3 at λ opt = 1028nm (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 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 FOM T at 0.913 is very close to the best 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 achieve over 5x enhancement of 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 pointsrandom 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 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 FOM T above the median, and the bottom inset shows the field-pattern with the highest FOM T . A comparison of all optimized field patterns at 1028nm and 625nm shows that the latter contains more higher-order Bessel-basis contributions.

Dependence on Illumination Wavelength
We further investigate the influence of λ opt by plotting optical force or torque per incident-field intensity I avg as a function of illumination wavelength. In Figs. 4(a) and 4(b), 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. 4(a)) and also enhance the force at the quadrupole mode while suppressing the dipole mode ( Fig. 4(b)).
In Figs. 5(a)-5(c), 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 λ opt value without suppressing the quadrupole resonance. The optimized total fields in Fig. 5(d) at 1028nm and 805nm both exhibit a 4π phase change around the circumference of the particle, resembling a quadrupole resonance. Total-field distributions of initial random field (left) and final optimized field (right) for the best optimized field configurations at three target wavelengths. Scalebar is 400nm.
In Fig. 5(c), 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. 5(d) 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 or 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.

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 where δ is a complex-valued random error bounded by the magnitude of the largest coefficient multiplied by the fractional weight 0 ≤ w ≤ 1. Fig. 6 shows the decrease in FOM T as a function of w. At λ 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, FOM T of the best optimized field at λ quad dropped Figure 6: Robustness of the optimized incident fields, quantified by the decrease in FOM T with respect to fractional random error added to the best (blue) and the median (red) optimized fields, at λ dip (right) and λ quad (left). The error bar represents the standard deviation for 100 samples. 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.

Conclusion
We present a numerical framework for computer optimization of structured illumination that maximizes optical force or 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.