Multiple criteria optimization of electrostatic electron lenses using multiobjective genetic algorithms

Designing electrostatic electron optical systems manually can become quite challenging when the number of electrodes increases. A large number of parameters are involved, such as the electrode thicknesses, the spacings between the electrodes, the aperture sizes of the electrodes, and the electrode voltages. These can all be varied to optimize the performance of the optical system, while simultaneously taking typical constraints into account, such as the maximum allowable field between electrodes to prevent discharges. It was demonstrated¹ that the design of such systems with a single objective function, i.e., a single requirement that needs to be fulfilled for optimized performance, can be automated using optimization techniques based on genetic algorithms.² A six-electrode electrostatic lens system was optimized to focus a primary electron (PE) beam to the smallest possible spot in a fixed image plane. The objective function in this case was the spot size, consisting of contributions from the spherical aberration and the chromatic aberration of the lens. To determine the lens properties and the objective function requires calculation of the lens field. This can be done most accurately using a finite element method, such as that offered by the commercial simulation package COMSOL. However, calculating the objective function using COMSOL is relatively slow (∼1 min for each system evaluation). This makes the optimization extremely time consuming. A faster but less accurate method is to calculate the axial potential of the lens using the second-order electrode method (SOEM)^3–6 (∼1 s for each system evaluation). In evolutionary algorithms (EA) such as genetic algorithm (GA), the objective function usually needs to be called numerous times in the optimization process, as many as hundreds or thousands of times. Using COMSOL then makes the optimization extremely time consuming.

A promising technique to reduce the computation time in such problems is the use of surrogate models,⁷ also known as surrogates. These are computationally inexpensive approximation models, employed to assist in the evaluation of computationally expensive functions to reduce the computational time of the problem. Studies on surrogate-assisted evolutionary algorithms such as GA began a decade ago. In the subsequent years, this method has been implemented in many different applications to optimize single and multiobjective optimization problems and found to be very efficient.^8–10 The surrogates (called the “low fidelity model”) can be implemented through different strategies combined with the accurate function calculation (known as the “high fidelity model”) to assist the optimization process.¹¹ 

In the previous study,¹ a surrogate-assisted GA was used by first running a GA optimization using field calculation by SOEM (low fidelity), resulting in a number of approximately good systems that were then inserted into a subsequent GA optimization using field calculation by COMSOL (high fidelity). With this approach, an optimized design was obtained in a reasonable time.

In this work, a similar optical system is considered but now with two objective functions. Challenging and labor-intensive examples of such electron optical designs can be found here,^12,13 but the example to be used in this work is a five-electrode lens that focuses a PE beam in a fixed sample plane and simultaneously projects the secondary electrons (SEs) emitted from the sample back through the lens onto a detector. This situation is typically encountered in electrostatic scanning electron microscopes with in-lens detection. To optimize this system requires a multiobjective optimization.

A recent work¹⁴ has been published on the design of an ion optical device with many possible free parameters, involving multiobjective functions (up to two). There, an adjoint variable method is implemented. In the optimization based on the adjoint method, the computational time would stay nearly constant with the increment of the free parameters, in contrast to the optimization using the evolutionary algorithms where it would increase dramatically. This gives an advantage over the evolutionary optimization technique for electron/ion design problems, where many free parameters are involved. However, the optimization based on the adjoint method is suitable for a few objective functions, while evolutionary optimization techniques such as GA are capable of optimizing multiple objective functions. Another difference is that the adjoint-based optimization is a local optimization technique (which can be used for global optimization problems by sampling different initial points), whereas the evolutionary algorithms are considered as global optimization techniques due to their metaheuristic characteristics that can automatically search more spaces through the objective function landscape and prohibit trapping in local optima.

In the design of electrostatic lens systems, such as in our case-study, the objective function landscape (even in the case of a single objective function) appears to have many local minima.¹⁵ In addition, multiple objective functions are usually aimed to be optimized, while the number of free parameters is not huge (in the order of tens). In such circumstances, EA/GA appears to be a better fit to the problem. Hence, a multiobjective genetic algorithm (MOGA)¹⁶ is proposed to be used here. A surrogate model, similar to that in the previous work,¹ is then applied in the problem to assist MOGA.

It is by no means obvious that the same two-step optimization strategy, combining SOEM and COMSOL, as was successfully used in the single objective function case, can be applied here too. The different energy and angular distributions of the PE and SE may cause the latter to follow trajectories considerably further away from the optical axis than the former. This means that the potential as approximated by the use of SOEM may not be sufficiently accurate to trace the SE toward the in-lens detector. Although a more accurate potential can be obtained using COMSOL, this will make the optimization really slow (one to two orders of magnitude slower, compared to optimization using SOEM), especially considering that a MOGA optimization scheme with two nonlinear objective functions needs thousands of runs.¹⁵ 

The objective of this work is to investigate whether the combination of an initial optimization using the less accurate SOEM and a subsequent optimization using COMSOL can also optimize systems with two objective functions in a reasonable time.

A schematic drawing of the electrostatic electron probe-forming objective lens used in this study is presented in Fig. 1. The lens, consisting of five planar electrodes, in combination with the sample positioned in a fixed image plane, is used to focus the PE beam, coming from the left, onto the sample plane. This geometry resembles an electrostatic scanning electron microscope design presented elsewhere,^17–19 and some of the dimensions were taken from there.

The SEs emitted from the sample by exposure to the PE beam are accelerated into the lens and directed toward the detector. The detector is a disk-shaped electrode (radius of 2.5 mm) with a central hole (radius of 0.5 mm) around the optical axis. The origin of the coordinate along the optical axis is taken at the surface of the first electrode, i.e., the surface closest to and facing the detector. The detector plane is positioned at −23 mm and the sample plane is at +17 mm.

During the optimization process, the parameters that define the geometry of the lens, such as the thickness of the electrodes (Ti), the aperture radius of the electrodes (Ri), and the gap between electrodes (Gi) are allowed to vary within certain boundaries. The allowed intervals are 1 mm < Ti (i = 1, …, 5) < 3 mm, 0.1 mm < Ri (i = 1, …, 5) < 2 mm, and 1 mm < Gi (i = 1, …, 4) < 3 mm. The voltage of the first electrode V1 is allowed to vary from 6 to 8 kV, whereas the other voltages Vi (i = 2, …, 5) can range from 700 to 10 kV. There are 19 free variables in total. The detector is kept at the same voltage as the first electrode V1. The voltage at the image plane is fixed to 600 V. Two constraints are set: a fixed image plane position (X-crossover = 17 mm) and a voltage breakdown condition (electric field < 15 kV/mm).

The first objective function is the probe size, and its calculation is described in detail in Ref. 1. In brief, the probe size contains contributions from two axial aberrations only (spherical and chromatic aberration) and is given by²⁰ 

where C_s and C_c are the spherical and chromatic aberration coefficients in the image plane, respectively; α is the half opening angle at the probe (here taken to be 10 mrad); U is the landing energy at the sample; and ΔU is the energy spread in the PE beam (here taken as 1 eV). In Eq. (1), the contributions from the geometric source image and diffraction are neglected. This is because, for very small probe currents, the contribution from the source image is negligible compared to those from the axial aberrations and the contribution from diffraction does not change during optimization, since the landing energy and the half opening angle of the PE beam at the probe are fixed. The spherical and chromatic aberration coefficients can be calculated from aberration integrals, which contain the axial potential, its derivatives, and a principle ray r_a(z), starting in the object on the optical axis with a unit angle (45°).²¹ Furthermore, the magnification has to be taken into account, which is obtained by tracing a second principle ray r_b(z) from the object plane, at unit height and zero angle, to the image plane. The position of the image plane is found where r_a(z) crosses the optical axis.

The second objective function is the detector collection efficiency of the SEs emitted from the sample. This requires tracing SE from the sample plane at z = +17 mm to the detector plane at z = −23 mm. The starting angles of the SE emerging from the sample are chosen from a uniform distribution of polar angles θ (the emission angle with respect to the sample surface normal) between 0.01 and 1.5 rad. The azimuthal angle is kept constant because of the rotational symmetry of the lens. However, in reality, the angular distribution of the SE yield is proportional to cos θ^22,23 and the SEs are emitted within a solid angle dΩ = 2π sin θ dθ. Therefore, the detector signal is weighed by a factor of sin θ cos θ ∼ sin 2θ. For simplicity, the energy of the emitted SE is taken from a uniform distribution between 1 and 10 eV.

The ray tracing and the determination of the lens properties require calculation of the lens field. The SOEM approach^1,3 makes use of the fact that the solution to the Laplace equation can be expressed in terms of the axial potential and its derivatives with respect to z. By ignoring terms higher than second order and using a cubic spline approximation to the axial potential, the latter can be obtained by solving a set of linear equations. The results are accurate within the paraxial approximation and prone to deviations further from the optical axis.^1,6 COMSOL is an accurate method to calculate the potential in the space of the lens system. However, as this method meshes the entire lens space to calculate the potential using the finite element method, it is associated with long computation times.¹ Both approaches will be used to calculate the objective functions of the example system.

The ray tracing is done using a MATLAB code. For the PE, which follows small angle trajectories, the equation of motion in the paraxial approximation is used, as in Ref. 1. However, for the SE, which follows larger angle trajectories, a more accurate real ray tracing is performed. The equations of motion are as follows:²⁴ 

Here, E_r and E_z are the radial and axial components of the field, respectively. The equations are solved numerically by taking very small time-steps of Δt, starting from (r₀, z₀), with initial velocity of v_0r and v_0z (calculated based on the initial energy of the SE).

In this section, the effect is studied, on both objective functions, of calculating the potential by SOEM, respectively, by COMSOL. Its effect on the probe size objective function can be judged from the two principle rays r_a(z) [i.e., the ray starting in the object on the optical axis with unit angle (45°)] and r_b(z) (i.e., the ray starting from the object plane at unit height and zero angle). Figure 2 shows, for a typical system within the range mentioned in Sec. II, both rays traced through the potential as calculated by SOEM, respectively, COMSOL. It is seen that the principle rays calculated by SOEM has a slight deviation but reasonably good overlap, with respect to the one calculated by COMSOL.

It is to be expected that the method of calculating the potential has a greater effect on the second objective function. For SE, starting at the sample plane (z = 17 mm) with angles ranging from 0.01 to 0.5 rad, and with an energy of 1 eV, the trajectories are calculated through fields calculated by SOEM and COMSOL to the detector at z = −23 mm. The results are shown in Fig. 3, where in Fig. 3(a) the SOEM-calculated trajectories are shown as dashed lines and the COMSOL-calculated trajectories as solid lines. Figure 3(b) shows the same rays on a larger vertical scale including the lens geometry. Although the SOEM rays clearly deviate from the COMSOL rays, they are not dramatically far off. Considering that the required computation time for SOEM-rays is around 1 s and for COMSOL rays around 60–70 s, it may be worthwhile in the optimization of lens systems to initially use SOEM for the potential calculation, thereby creating some approximately good systems, and feed those to a more accurate COMSOL-based multiobjective function genetic algorithm (MOGA) optimization. This is the subject of Sec. V.

The GA, categorized as one of the evolutionary algorithms, mimics natural evolution (inspired by Darwin’s theory). It starts with a randomly generated initial population that includes a set of systems, denoted by P₁(x₁, … , x_n). In natural evolution, “x_i” represents the chromosomes of different members. In electron lens system optimization, “x_i” represents the electron lens systems. Hence, P₁ includes a set of initially randomly generated electron lens systems, defined based on the lens system variables. In nature, across time, the initial population is gradually improved toward members that are better matched with their environmental conditions. In GA, analogously, the initial population evolves toward a new set of systems [P_i+1(x₁, … , x_n)] which better satisfy the conditions of the problem in each so-called “generation.” The new population is mainly created by the operators “selection,” “crossover,” and “mutation.” Conditions that are determined to be optimized are formulated by a so-called objective function.

The multiobjective GA (MOGA)¹⁶ differs from a classical GA in how the objective function value is assigned and ranked to each member in the population. The initial population is created randomly, similar to that for the GA. The next generations are derived using the nondominated ranking.²⁵ Each member gets a nondominated rank corresponding to its relative objective function. The members of the population are then evaluated, ranked, sorted, and selected by the MOGA based on the nondominated classification. A set of pareto fronts²⁵ (the nondominated solutions for which none of the objective functions can be improved in value without degrading the values of other objective functions) are found and presented to the user at the end of the optimization process. The user has then the choice to select the system that best suits their problem among them. The remaining features of the algorithm, such as “crossover” and “mutation,” are the same as those in a classical GA.

The proposed algorithms here are run on a PC with an Intel (R) Xeon (R) W-2123 CPU @3.60 GHz and 32 GB of RAM.

It will first be analyzed how MOGA optimization of the example lens system performs when SOEM is used for the potential calculation (hereafter called MOGA + SOEM). The parameters used for the MOGA in this case study are population size = 100, maximum generations = 1000, and crossover fraction = 0.5. The optimization starts with an initial population of 100 systems randomly created by MOGA. Multiple runs are performed from which two example runs are presented in Figs. 4(a) and 4(b). The figures show for a typical run the SE detection efficiency versus the probe size. In a relatively short time of 20 h, MOGA evaluated 100 000 systems from which 30 000–40 000 systems passed the constraints. Hereafter, these systems are referred to as good systems.

In Figs. 4(a) and 4(b), which do not include all systems but only the good systems, a trend towards decreasing probe size and increasing detection efficiency is seen. The pareto front²⁵ points for the two runs are illustrated as the dark blue and cyan points. By MOGA optimization of the voltages and the geometry of the lens, the initial probe size of 60 nm is reduced to only 9 nm with a detection efficiency ranging from 0% to 49%. The system with a probe size of 18 nm and a detection efficiency of 49% is indicated by “L” in Fig. 4(a), the system with 9 nm probe size and 14% detection efficiency by “K” in Fig. 4(b). It takes on average ∼200 min to reach the pareto front systems in this optimization.

Figure 4(c) shows Figs. 4(a) and 4(b) in one graph. It is clearly seen that for different runs the pareto fronts vary considerably. This reveals the complexity of the (highly nonlinear) objective functions and the existence of multiple local minima. The large number of free parameters in the design problem can easily cause the optimization to reach a different local optimum (pareto front) in different runs.

Next, a MOGA optimization is performed where for the calculation of the objective functions the necessary lens potential is calculated using COMSOL. Hereafter, we will call this method MOGA + COMSOL. All systems in the initial population are created randomly by MOGA itself. Since the computation time for the objective function calculation by COMSOL is dramatically longer, the population size is taken here to be 50 instead of 100, and the maximum generations as 200 instead of 1000. The crossover fraction is kept the same as before, i.e., 0.5.

In total, 10 000 systems have been evaluated and the entire optimization took about 7000 min (∼5 days). Of these 10 000 systems, 3000–4000 systems passed the constraints and are thus labeled as good systems. Figures 5(a) and 5(b) show the detection efficiency versus probe size for the good systems of two separate runs. The pareto front points of these two runs are shown as the black and orange points. Figure 5(c) shows the overlap of the plots in Figs. 5(a) and 5(b).

Table I summarizes the optical parameters for three different systems in the pareto-front of these MOGA + COMSOL optimizations. As can be seen from Table I, it takes at least ∼900 min of optimization time to reach the pareto-front system with a probe size of 17.5 nm with an SE detection efficiency DE_SE = 57% [system O in Fig. 5(c)] or ∼2000 min to reach a system with spot size D_spot = 13.0 nm and DE_SE = 52 [system N in Fig. 5(c)]. The smallest spot size reached in this optimization was D_spot = 12.3 nm with DE_SE = 48% nm [system M in Fig. 5(c)], which took ∼4000 min. It can be concluded that the computational time of MOGA + COMSOL to reach a satisfactory result is extremely long. In the next subsection, the use of good-systems from MOGA + SOEM as an input for MOGA + COMSOL is considered as a means to reduce the optimization time.

In an attempt to reduce the computational times involved in MOGA + COMSOL optimization, the initial input systems with preoptimized geometries and voltages obtained by MOGA + SOEM were fed into MOGA + COMSOL. Hereafter, this will be called MOGA + SOEM + COMSOL. 10 000 systems with the same optimization parameters as used in MOGA + COMSOL (population size = 50, maximum generations = 200, and crossover fraction = 0.5) were evaluated. 20 systems are taken from pareto front systems found by MOGA + SOEM and fed into the initial population, together with 30 systems created randomly by MOGA. The MOGA starts with these 50 initial systems while calculating the objective functions based on potential calculation by COMSOL. This is schematically shown as a flow chart in Fig. 6.

The optimization results for two runs are plotted as SE detection efficiency versus probe size in Figs. 7(a) and 7(b). The pareto front systems are indicated as well by the red and green stars in Figs. 7(a) and 7(b), respectively. The results from both runs are overlapped in Fig. 7(c). Table II presents the optical parameters of three nondominated optimized systems taken as representatives from the pareto front systems. It is seen that after only ∼100 min, the optimization found an optimized system having a probe size of 17.5 nm and a very high detection efficiency of 85% [system R in Fig. 7(c)], or a system with a probe size of 12.2 nm and a detection efficiency of 81% [system Q in Fig. 7(c)], found after ∼300 min. After ∼2500 min, this optimization reached a system with the very small spot size of 8.6 nm and detection efficiency of 73% [system P in Fig. 7(c)].

It is noted that the time consumed in MOGA + SOEM to generate the initial data (∼200 min) should be added to the computational time of this optimization.

Comparing the data in Tables II and I, the MOGA + SOEM + COMSOL optimization achieved much better systems in a considerably shorter time than the MOGA + COMSOL optimization. For instance, reaching a system with a spot size of 17.5 nm only takes 100 min in the former case and 900 min in the latter, even with a dramatically lower SE detection efficiency. It should be noted that the minimum spot size of 8.6 nm with a detection efficiency of 85%, reached in the time frame of ∼300–2500 min, could not be reached at all by the MOGA + COMSOL optimization, even after the evaluation of 10 000 system (∼7000 min).

In Fig. 8, the results from Figs. 5(c) and 7(c) are compared in one graph, i.e., a comparison between optimized systems obtained by MOGA + COMSOL and MOGA + SOEM + COMSOL, respectively, including the pareto front systems. It is noted that the runs in the latter case show much more overlap than what was seen in the former case. Furthermore, the systems optimized by MOGA + SOEM + COMSOL clearly outperform the ones optimized by MOGA + COMSOL. Hence, it is concluded that, despite the deviating field calculation by SOEM, it is still very effective to use this technique in the optimization, to find preliminary, approximately good-systems, to be fed into another optimization, which is performed based on the accurate field calculation by COMSOL, and to boost the optimization.

The evolution of some lens systems during the optimization by MOGA + SOEM + COMSOL is summarized in Table III. A clear change is observed in the geometries and voltages of the systems while the optimization is progressing until it finally reaches a set of optimized systems. It should be noted that those SE hitting the electrodes do not contribute to the detection efficiency. The PEs are started at z = −23 mm, with angles ranging between 0.1 and 0.8 rad, and traced by paraxial ray-tracing. The horizontal axis runs from −10 to +17 mm for a better visualization. SE trajectories are produced by real ray-tracing. In the figures of Table III, these are drawn for visualization purposes only.

To validate the final optimization results, an optimized system found by MOGA + SOEM + COMSOL (system A4 in Table III) is simulated using the EOD software package from SPOC.²⁶ Figure 9 shows the EOD simulation result for the geometry and the trajectories of both PE and SE in this system. The optical parameters of system A4, as obtained from MOGA + SOEM + COMSOL and EOD, are compared in Table IV. The presented results show that the optimization result agrees very well with the EOD simulation.

The design of a pure electrostatic lens system with five electrodes and two objective functions was optimized using a surrogate-assisted multiobjective genetic algorithm (MOGA) approach. The system is required to focus a PE beam into a fixed image plane with minimum probe size and attract SE emitted from the image plane back into the lens on to an SE detector with the highest possible detection efficiency. The main challenge of such optimization is the long computation time necessary for an accurate COMSOL-based field calculation. The much faster, but less-accurate, field calculation technique SOEM combined with MOGA does not provide the best optimized systems, but reasonably good ones that can be used as input for a COMSOL-based MOGA optimization. This strategy has been demonstrated to lead to optimized systems in a reasonably short time of a few hours compared to many days when using COMSOL-based MOGA only. The optimization process is fully automated and may help electron optics designers to quicker find solutions for complex electrostatic systems.

The data that support the findings of this study are available within the article.