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ⁱⁿ lies in the x-y-plane; see Fig. 1.

 

Fig. 1 A dielectric lens D_L is illuminated by a plane wave Eⁱⁿ. The goal is to find the shape of D_L so that the focused light in D_F is maximized.

Download Full Size | PDF

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]

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)

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 ℝ². It is well known that T_V is a diffeomorphism if ‖V‖_{C¹(ℝ²;ℝ²)} < 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

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

 

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.

Download Full Size | PDF

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]

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

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 ${\text{V}}_N^{\text{update}}$ 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.

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

 

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.

Download Full Size | PDF

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² = 49 (b), 17² = 289 (c), and 27² = 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].

 

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.

Download Full Size | PDF

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.

 

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.

Download Full Size | PDF

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.