## Abstract

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

## 1. 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.

## 2. Model problem

Let *D*_{L} be the cross section in the x-y-plane of a simply connected, homogeneous, non-magnetic, 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
${\mathbf{E}}^{\text{in}}(\mathbf{x}):={E}_{z}^{\text{in}}{\mathbf{e}}_{z}\text{exp}(i\mathbf{k}\cdot \mathbf{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.

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]

*E*can be computed employing the finite element method [7]. Firstly, we introduce a bounded domain Ω that encloses

_{z}*D*

_{L}. Then we state the variational formulation of Eq. (1) on Ω [7]. That is,

*E*must satisfy

_{z}*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})

*D*

_{L}is optimal if it solves the following PDE-constrained shape optimization problem

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 ℝ^{2}. It is well known that *T*_{V} is a diffeomorphism if ‖V‖_{C1(ℝ2;ℝ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

**M**

_{V}is given by ${\mathbf{M}}_{\text{V}}=\left(\text{det}\mathbf{D}{T}_{\text{V}}\right)\mathbf{D}{T}_{\text{V}}^{-1}\mathbf{D}{T}_{\text{V}}^{-T}$, and

**D**

*T*

_{V}=

**I**+

**D**V, where

**I**is the identity matrix and

**D**V 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.

## 3. 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

*denotes the scalar*

_{i}*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

To compute the approximate solution ${E}_{z}^{h}$ 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]

*p*is the finite element solution of the adjoint problem

^{h}*χ*

_{DF}denotes the characteristic function of

*D*

_{F}.

Optimization is carried out iteratively by computing a descent direction ${\text{V}}_{N}^{\text{update}}$ as the solution of the linear variational problem

_{H1}we denote the usual

*H*

^{1}-inner product [7], and adding the update $\delta {\text{V}}_{N}^{\text{update}}$ to the map

*T*

_{VN}= I + V

*. The optimization step length*

_{N}*δ*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*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 ${\text{V}}_{N}^{\text{update}}$ it is necessary to evaluate

_{h}*d*J(V, ·) Eq. (8) for all basis vector fields B

*(*

_{i}**x**)

**e**

*and B*

_{x}*(*

_{i}**x**)

**e**

*. 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].*

_{y}Without Formula Eq. (8) at hand it would be virtually impossible to employ a steepest descent direction algorithm. Approximations of the gradient *d*J for 2*N* design parameters (see Eq. (7)) with finite differences would require at least 2*N* additional evaluations of the functional J (and thus to solve Eq. (6b) 2*N* additional times), whereas Formula Eq. (8) requires only the solution of Eq. (9). Note that 2*N* 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.

## 4. 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 optimize the shape employing 49 (b), 289 (c), and 729 (d) multivariate B-spline. Although the difference between the optimized and the initial shapes are hardly visible, we observe a high sensitivity of the field distribution. The target functional J increases to 162%, 203%, and 186% of the initial value, respectively. The improvement employing 789 B-splines is slightly less than the one for 289 B-splines. Again, this is a consequence of the ill-posedness of the optimization problem: there are many local minima close to the optimum, and steepest descent methods are not global minimizers [13].

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 non-smooth 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 B-splines 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.

## 5. 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.

## Acknowledgments

The work of A. Paganini and S. Sargheini was partly supported by ETH Grant CH1-02 11-1.

## 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). [CrossRef] [PubMed]

**2. **M.-S. Kim, T. Scharf, S. Mühlig, C. Rockstuhl, and H. P. Herzig, “Engineering photonic nanojets,” Opt. Express **19**, 10206–10220 (2011). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**4. **C. Hafner, “Boundary methods for optical nano structures,” physica status solidi (b) **244**, 3435–3447 (2007). [CrossRef]

**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). [CrossRef]

**7. **D. Braess, *Finite elements. Theory, Fast Solvers, and Applications in Elasticity Theory* (Cambridge University, 2007). [CrossRef]

**8. **J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**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). [CrossRef]