Shape optimization of microlenses

Microlenses are highly attractive for optical applications such as super resolution and photonic nanojets, but their design is more demanding than the one of larger lenses because resonance effects play an important role, which forces one to perform a full wave analysis. Although mostly spherical microlenses were studied in the past, they may have various shapes and their optimization is highly demanding, especially, when the shape is described with many parameters. We first outline a very powerful mathematical tool: shape optimization based on shape gradient computations. This procedure may be applied with much less numerical cost than traditional optimizers, especially when the number of parameters describing the shape goes to infinity. In order to demonstrate the concept, we optimize microlenses using shape optimization starting from more or less reasonable elliptical and semi-circular shapes. We show that strong increases of the performance of the lenses may be obtained for any reasonable value of the refraction index. © 2015 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.0180) Microscopy; (290.0290) Scattering; (080.4225) Nonspherical lens design. References and links 1. Z. Chen, A. Taflove, and V. Backman, “Photonic nanojet enhancement of backscattering of light by nanoparticles: a potential novel visible-light ultramicroscopy technique,” Opt. Express 12, 1214–1220 (2004). 2. M.-S. Kim, T. Scharf, S. Mühlig, C. Rockstuhl, and H. P. Herzig, “Engineering photonic nanojets,” Opt. Express 19, 10206–10220 (2011). 3. Y. Shen, L. V. Wang, and J.-T. Shen, “Ultralong photonic nanojet formed by a two-layer dielectric microsphere,” Opt. Lett. 39, 4120–4123 (2014). 4. C. Hafner, “Boundary methods for optical nano structures,” physica status solidi (b) 244, 3435–3447 (2007). 5. R. Hiptmair, A. Paganini, and S. Sargheini, “Comparison of approximate shape gradients,” BIT Numerical Mathematics, 1–27 (2014). 6. D. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory (Springer, 2013). 7. D. Braess, Finite elements. Theory, fast solvers, and applications in elasticity theory (Cambridge University, 2007). 8. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). 9. G. Allaire, Conception optimale de structures (Springer-Verlag, 2007). 10. K. Höllig and J. Hörner, Approximation and modeling with B-splines (Society for Industrial and Applied Mathematics, 2013). 11. M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich, Optimization with PDE constraints (Springer, 2009). 12. R. Hiptmair and A. Paganini, “Shape optimization by pursuing diffeomorphisms,” SAM-Report 2014-27, ETHZ (2015). 13. J. Nocedal and S. J. Wright, Numerical optimization (Springer, 2006). 14. K. Eppler and H. Harbrecht, “Coupling of FEM and BEM in shape optimization,” Numer. Math. 104, 47–68 (2006).


Introduction
It is well known that the classical diffraction limit of optical lenses restricts the resolution of optical microscopes to approximately half of the wavelength.One way to increase the resolution is to work within media (e.g., oil) with an index of refraction higher than the one of air, i.e., with a wavelength smaller than the one of air.Since lenses always have an index of refraction higher than the one of air, sufficiently thick lenses will have a focus point inside that is smaller than half the wavelength in air.By tuning the index of refraction of a microlens of spherical shape (3D) or of circular-cylinder shape (2D), one may move the focus on or close to the surface of the lens.By doing this, one may overcome the diffraction limit in some sense.Ultramicroscopy based on this concept has been proposed in [1].
When the focus of a microlens is on or near its surface, one may also observe a very much elongated shape of the focus area, which is called nanojet.Nanojets have been shown both by simulations and experiments [2].
The length of a nanojet may be strongly increased, for example, by considering layered spheres [3].It should be mentioned that these ultra-long nanojets have a waist of more than half wavelength, but they have various promising applications as mentioned in [3].
The main advantage of studying microlenses of spherical or of circular-cylinder shape is simplicity, i.e., only two parameters (radius and index of refraction of the lens) have to be studied for an incident plane wave of a certain wavelength.Furthermore, Mie scattering theory may be used for obtaining solutions with high accuracy.However, for engineering the shape of the nanojet, more complex models with more tuning parameters are needed.Furthermore, it is important to note that appropriate low-loss dielectrics for manufacturing microlenses are only available with certain values of the index of refraction, i.e., this index should not be considered as a parameter that may be optimized.
The standard way to tune the optical features of lenses is shape optimization.An obvious geometry with still very few optimization parameters is a 2D lens of cylindrical shape with elliptical cross section [4] or a 3D lens with the shape of an ellipsoid.In [4], it was also shown that lenses with another topology (e.g., toroidal) may cause nanojets that are located considerably away from the surface.
By increasing the number of parameters describing the microlens, one can increase the search space for finding nanojets with desired shape and one may hope for better solutions than what has been found so far.This leads to very demanding optimization problems because the computation time of conventional numerical optimizers increases drastically with the number of parameters.In this paper, we demonstrate that this problem may be overcome by using advanced mathematical techniques, so-called shape gradients [5].We outline the concept of shape optimization.Based on shape gradients, this allows us to start with a lens of arbitrary shape and index of refraction and to deform it until a local optimum is reached.We can handle an unlimited number of shape parameters and obtain various promising structures without high computational cost.However, it goes without saying that best results are achieved with a reasonable initial guess.
In order to demonstrate the procedure, we consider a rather simple 2D test problem, starting from elliptic semicircular lenses.We optimize thier shape in such a way that the field inside a specified rectangular box is maximized.This procedure can easily be applied to lenses of axisymmetric shape and for various optimization goals, for example, strong gradients of the intensity near the nanojet, which is interesting, e.g., for optical trapping.Since we utilize a standard finite element method (FEM), a full 3D particle shape could also be considered.
Let D L be the cross section in the x-y-plane of a simply connected, homogeneous, nonmagnetic, cylindrical lens that is infinitely long in the z-direction and has a refractive index n L .The lens is surrounded by non-magnetic material with refractive index n A = 1.It is illuminated by a plane wave E in (x) := E in z e z exp(ik • x), where e z denotes the unit vector in the z-direction.The wave vector k of the incident field E in lies in the x-y-plane; see Fig. 1.
Fig. 1.A dielectric lens D L is illuminated by a plane wave E in .The goal is to find the shape of D L so that the focused light in D F is maximized.
Because of translation invariance of electric and magnetic fields in the z-direction, and tangential continuity of the magnetic field on ∂ D L , the z-component of the total electric field E z satisfies the partial differential equation (PDE) [6] Additionally, the scattered field Approximations of E z can be computed employing the finite element method [7].Firstly, we introduce a bounded domain Ω that encloses D L .Then we state the variational formulation of Eq. ( 1) on Ω [7].That is, E z must satisfy where H 1 (Ω) denotes the Sobolev space of square integrable functions with square integrable weak derivatives [7].Finally, we replace the right hand side of Eq. ( 3) by realizing the radiation condition Eq. ( 2) via a perfectly matched layer (PML) [8].
We want to find the shape of the lens D L so that the focused light in the focus or nanojet area D F is maximized.To model this target optical property, we introduce the shape functional (energy content in D F ) Then, we say that the shape of the lens D L is optimal if it solves the following PDE-constrained shape optimization problem argmax To construct the set of admissible shapes U ad (D L,0 ), we collect all domains that can be obtained by letting a diffeomorphism (a bijective continuously differentiable map with continuosly differentiable inverse) T act on an initial guess D L,0 for which Eq. ( 1) is well-defined.In particular, we consider maps of the form T V (x) := x + V(x), where V is a vector field on R 2 .It is well known that T V is a diffeomorphism if V C 1 (R 2 ;R 2 ) < 1 [9].Note that all admissible shapes share the same topology.
Since T V is a diffeomorphism, we can employ transformation techniques (change of variables) to replace the dependence of Eq. ( 5) on D L with a dependence on V. We obtain the equivalent optimal control problem argmax where the matrix M V is given by and DT V = I + DV, where I is the identity matrix and DV is the Jacobian of V.With this approach, it is sufficient to construct an initial grid on D L,0 to simulate any shape in U ad (D L,0 ), avoiding the need to generate a new FEM mesh for every design.

Discretization and optimization
To discretize the control V we employ multivariate B-splines of degree 3 [10] generated on a regular grid that covers D L and does not intersect D F ; see Fig. 2.More precisely, we consider vector fields that can be written as where B i denotes the scalar i-th multivariate B-spline of degree 3. A multivariate B-spline is the tensor product of univariate B-splines, which are piecewise polynomials with compact support.For example, the formula of a univariate cubic B-spline on the knot sequence (0, 1, 2, 3, 4] reads This combination of features is excellent to access local sensitivity information and allows for an efficient implementation. To compute the approximate solution E h z of the state constraint Eq. (6b) we employ linear Lagrangian finite elements on quasi-uniform triangular meshes [7].
Since both the control V and the state variable E z belong to infinitely dimensional spaces, the discrete counterpart of Eq. ( 6) is inherently a high dimensional optimization problem.To solve it we employ a steepest descent direction algorithm.
The derivative in the direction W of the target functional J evaluated in V can be computed following the standard differentiation techniques for PDE-constrained functionals [11].By the nature of the model problem under consideration it is meaningful to assume that the support of the vector fields V and W does not intersect the region of interest D F .Then, the derivative of J in the direction W reads [12]  where The scalar function p h is the finite element solution of the adjoint problem Sommerfeld radiation condition Eq. ( 2), where χ D F denotes the characteristic function of D F .Optimization is carried out iteratively by computing a descent direction V update N as the solution of the linear variational problem where by (•, •) H 1 we denote the usual H 1 -inner product [7], and adding the update δ V update N to the map T V N = I + V N .The optimization step length δ is computed according to [13].
An optimization step comprises the computation of the finite element solution of Eq. (6b) and Eq. ( 9), and of the descent direction in Eq. (10).To compute u h and p h one needs to assemble a (sparse) stiffness matrix (the same can be used for both) and to solve the related linear system.The latter task can be performed efficiently exploiting sparsity, whilst the computational cost of matrix assembly is directly proportional to the number of triangles of the mesh, and is only slightly larger than the assembly of the stiffness matrix related to Eq. (1).To compute the descent direction V update N it is necessary to evaluate dJ(V, •) Eq. ( 8) for all basis vector fields B i (x)e x and B i (x)e y .Although with a larger constant than the one for the matrix assembly, its computational cost is directly proportional to the number of triangles of the mesh because B-splines have compact support.A detailed explanation of an efficient implementation of the algorithm in MATLAB is provided in [12].
Without Formula Eq. ( 8) at hand it would be virtually impossible to employ a steepest descent direction algorithm.Approximations of the gradient dJ for 2N design parameters (see Eq. ( 7)) with finite differences would require at least 2N additional evaluations of the functional J (and thus to solve Eq. (6b) 2N additional times), whereas Formula Eq. (8) requires only the solution of Eq. ( 9).Note that 2N additional evaluations of J only lead to a first order approximation of the gradient, which would be considerably less accurate than the proposed method.Finally, due to the high number of design parameters, performing optimization with a genetic algorithm is not computationally affordable because it would require too many evaluations of the functional J.

Numerical experiments
We first start with an initial guess D L,0 that is an ellipse with semi-minor and semi-major axes of length 4λ and 5λ , respectively.We illuminate D L,0 from the left as shown in Fig. 1 and locate the region of interest D F , a rectangle of sides λ /2 and 2λ , on the backside of D L,0 .The boundary of the grid used to generate B-splines is a square of sidelength 11.2λ as in Fig. 3.The mesh employed for finite element computations has 211313 nodes, 421184 triangles, and width h = 0.083λ .In the first numerical experiment we choose the refractive index n = 1.5.Before optimization, the focus lies beyond the region of interest; see Fig. 3(a).We optimize the shape employing 7 2 = 49 (b), 17 2 = 289 (c), and 27 2 = 729 (d) multivariate B-splines.The target functional J increases to 315%, 355%, and 359% of the initial value, respectively.Essentially, the optimized shapes are thicker than the original ellipse.This obviously shifts the focus close to the lens surface.The three retrieved shapes differ slightly, but this is a consequence of the optimization problem being ill-posed [14].In particular, we observe that the front of the lens evolves into a more or less sinusoidal line, when sufficiently many B-splines are used.
In the second numerical experiment we choose the refractive index n = 2.7.In this case, the initial shape is already almost optimal; see Fig. 4(a).As in the first numerical experiment, we In the third numerical experiment we choose the refractive index n = 3.5.Before optimization, the focus lies inside the lens; see Fig. 5(a).As in the first numerical experiment, we optimize the shape employing 49 (b), 289 (c), and 729 (d) multivariate B-spline.This time, optimized shapes are thinner in order to shift the focus outside the lens.In this case, the high sensitivity of the field distribution with respect to the shape boundary influences the improvements of the target functional J, which increases to 451%, 461%, and 570% of the initial value, respectively.
Finally, we perform a fourth numerical experiment to show that we can easily deal with nonsmooth shapes, too.The initial guess D L,0 is a half-circle of radius 2λ with refractive index n = 2; see Fig. 5(e).The boundary of the grid used to generate B-splines is a square of side 5λ as in Fig. 5(e).The mesh employed for finite element computations has 211937 nodes, 422432 triangles, and width h = 0.0817λ .The optimized shape obtained employing 729 Bsplines is displayed in Fig. 5(f).The target functional J increases to 217% of the initial value.Additionally, we observe that the retrieved shape reflects less.We repeat the experiment for n = 1.5 and n = 3.The optimized shapes are displayed in Fig. 5 ((g),(h), respectively).The target functional J increases to 183% and 482% of the initial value, respectively.For n = 1.5 the thickness of the lens in the region in front of D F is increased, whilst for n = 3 optimality is achieved by corrugating the surface so that transmission increases (see Fig. 5(i)), i. e., the optimization runs towards a mixture with a Fresnel-type lens.

Conclusion and outlook
We have outlined and applied shape optimization.In order to demonstrate the concept, numerical models of microlenses with different indices of refraction, starting from elliptical and semi-circular shapes, were optimized and various shapes were obtained.The goal of these optimizations was to maximize the energy in a rectangular output box behind the lens.It should be pointed out that one may apply the procedure to problems with other optimization goals, e.g., maximum intensity gradients in the output box, which would be interesting for optical tweezers; other shapes of the output box; more than one output box; etc.Interestingly, the reflection of the microlenses was also reduced although this goal was not explicitly formulated.The reason for this is that increasing reflection will decrease the light intensity in the output box.For a proof of concept, only 2D lenses were considered, but the procedure can easily be extended to more realistic axisymmetric lenses or even, with some bigger numerical effort, to full three-dimensional lenses.Note that shape optimization is not limited to the optimization of lenses.It could also be applied to various other applications, among which, for instance, plasmonic antennas or optical tweezers.

Fig. 2 .
Fig. 2. Left: Grid used to generate multivariate B-splines of degree 3. Right: Multivariate B-splines of degree 3. Its support comprises 4 × 4 grid cells.The B-spline is polynomial in each cell.

Fig. 3 .
Fig. 3. First numerical experiment: Absolute value of E z before (a) and after optimization with 49 (b), 289 (c), and 729 (d) multivariate B-splines.The optimized shapes are thicker in order to shift the focus close to the lens surface.

Fig. 4 .
Fig. 4. Second numerical experiment: Absolute value of E z before (a) and after (b,c,d) optimization.We observe a high sensitivity of the field distribution.

Fig. 5 .
Fig. 5. Third numerical experiment: Absolute value of E z before (a) and after (b,c,d) optimization.The optimized shapes are thinner in order to shift the focus outside the lens.Fourth numerical experiment:Absolute value of E z before (e) and after (f) optimization for n = 2, and upper half of optimized and initial lens for n = 1.5 (g) and n = 3 (h).For the latter case, we plot |E z | along the x-axis before (red dotted line) and after (blue line) optimization (i).We observe that the optimum is achieved by drastically increasing transmission.