A new approach for fast field calculation in electrostatic electron lens design and optimization

In electron optics, calculation of the electric field plays a major role in all computations and simulations. Accurate field calculation methods such as the finite element method (FEM), boundary element method and finite difference method, have been used for years. However, such methods are computationally very expensive and make the computer simulation challenging or even infeasible when trying to apply automated design of electrostatic lens systems with many free parameters. Hence, for years, electron optics scientists have been searching for a fast and accurate method of field calculation to tackle the aforementioned problem in the design and optimization of electrostatic electron lens systems. This paper presents a novel method for fast electric field calculation in electrostatic electron lens systems with reasonably high accuracy to enable the electron-optical designers to design and optimize an electrostatic lens system with many free parameters in a reasonably short time. The essence of the method is to express the off-axis potential in an axially symmetrical coordinate system in terms of derivatives of the axial potential up to the fourth order, and equate this to the potential of the electrode at that axial position. Doing this for a limited number of axial positions, we get a set of equations that can be solved to obtain the axial potential, necessary for calculating the lens properties. We name this method the fourth-order electrode method because we take the axial derivatives up to the fourth order. To solve the equations, a quintic spline approximation of the axial potential is calculated by solving three sets of linear equations simultaneously. The sets of equations are extracted from the Laplace equation and the fundamental equations that describe a quintic spline. The accuracy and speed of this method is compared with other field calculation methods, such as the FEM and second order electrode method (SOEM). The new field calculation method is implemented in design/optimization of electrostatic lens systems by using a genetic algorithm based optimization program for electrostatic lens systems developed by the authors. The effectiveness of this new field calculation method in optimizing optical parameters of electrostatic lens systems is compared with FEM and SOEM and the results are presented. It should be noted that the formulation is derived for general axis symmetrical electrostatic electron lens systems, however the examples shown in this paper are with cylindrical electrodes due to the simplicity of the implementation in the software.

In electron optics, the quality of an electrostatic lens system is specified by its electron-optical properties, which are determined from its electric field 1-3 .Calculation of the electric field for electrostatic lens systems is hence the most significant step needed in evaluation of the lens system.Therefore, the speed and accuracy of calculating the electric field (or electric potential) plays a major role in electron-optics computation, design and optimization.The electric field can be calculated using numerical methods such as FEM 4 , Boundary Element Method (BEM) 5  and Finite Difference Method (FDM) 6 .Although these methods result in very accurate values for the electric field (which consequently results in accurate determination of optical parameters), all these methods are time consuming.For instance, calculation of the potential distribution using FEM and performing ray tracing takes ̴ 1 min for each system 7 .Note: All computer simulations in this study were conducted on our system, a PC equipped with an Intel® Xeon® W-2123 CPU @3.60 GHz and 32 GB of RAM, using the MATLAB programming language.The reported computation times are estimated based on the results obtained on the mentioned PC and www.nature.com/scientificreports/Additionally, a comparison with the previously introduced SOEM based optimization method (SOEM + FEM + GA 7,12 ) is essential to investigate the superiority of the newly introduced method based on FOEM in optimizing electrostatic lenses.
This paper is organised as follows.In "Mathematical derivation of the axial potential for electrostatic electron lens systems using FOEM" section, the mathematical derivation for the axial potential calculation using FOEM is explained and presented.To derive the axial potential formulation, first the fundamental equations of the quintic spline 26 for a dataset with unequal gaps between the data points are derived in "Fundamental equations of the quintic spline".Next, in "Solving the Laplace equations up to the fourth order derivatives" section, using these equations in combination with the Laplace equation and proper boundary conditions, a system of linear equations is formed, from which the axial potential and its first and second derivatives are calculated.The accuracy and speed of the field and optical parameters calculation of this new method in comparison with FEM and SOEM is presented in "Accuracy of FOEM in field and optical parameters calculation" section for different typical electrostatic lens systems including 3, 4 and 5 lens electrodes.In "Optimization of electrostatic electron lens systems based on FOEM and a comparison with SOEM-and FEM-based optimization" section, the newly introduced field calculation method (FOEM) is implemented in a genetic algorithm-based optimization program (FOEM + GA).
Subsequently, in "The calculation of optical parameters in a single run of GA Optimization based on FOEM, versus SOEM and FEM" section, the improvement trends and accuracy calculations of optical parameters for the optimized systems in a single run of GA-aided optimization based on FOEM are investigated and compared with situations where the optimized systems are calculated using SOEM and FEM.This is done to obtain more statistically reliable results for the accuracy comparisons of optical parameters calculated using different methods of FOEM, FEM, and SOEM.Additionally, we illustrate and compare the trends of improvements in finding better lens systems in the mentioned single run of GA optimization based on FOEM with situations where the optimized systems are calculated using SOEM and FEM.
In "The efficiency of GA optimization based on FOEM, FEM and SOEM" section, a comparative study on the efficiency of various GA-aided optimization approaches based on FOEM, FEM, and SOEM is presented.Building on this study, to assess the efficiency of FOEM + GA, the results of FOEM + GA are fed into another optimization based on FEM (FOEM + FEM + GA), and the outcomes are presented at the end in "The efficiency of GA optimization based on FOEM, FEM and SOEM" section.Conclusions are drawn in "Conclusion" section.
Note: The formula derived and applicable for electrostatic electron lens systems in both cylindrical and noncylindrical configurations.However, the examples showcased in this paper utilize cylindrical electrodes.This choice is driven by the software's ease of implementation.

Mathematical derivation of the axial potential for electrostatic electron lens systems using FOEM
In this section, the aim is to derive the axial potential formulation for rotationally electrostatic lens systems with multi-electrodes.An illustrative example of such lens systems is provided schematically in both 3D (Fig. 1a) and 2D (Fig. 1b).Where, T j , R j , and V EL j refer to the thicknesses, radii, and voltages at each electrode, respectively ( j = 1, 2, . . ., tot , where tot is the total number of electrodes).While G j represents the gaps between two consecutive electrodes (( j = 1, 2, . . ., tot − 1).
It can be demonstrated that for such geometries with axial symmetry, the spatial potential distribution can be obtained from the axial potential and its derivatives.The solution of the Laplace equation for a rotationally symmetrical geometry in terms of the axial potential and its derivatives is given by 27   in which r is the radial distance to the axis, z is the horizontal distance from the entrance of the lens system, ϕ(r, z) is the potential in space at coordinate (r, z) , ϕ(0, z) is the axial potential at coordinate (0, z) , and ϕ (2i) (0, z) is the www.nature.com/scientificreports/(2i) th derivative of the axial potential at coordinate (0, z) .When the symmetry axis is discretized into N points, and the terms are truncated after the fourth order derivative of the axial potential, Eq. (1) can be expressed as: Considering Eq. ( 2), it is evident that the potential distribution along the axis ( ϕ i ) is linked to voltages at the surface of the electrodes ( V i (r i )) , where the values of V i (r i )) are known at each electrode and correspond to the electrode voltages assigned to that specific electrode ( V ELj , j = 1, 2, . . ., tot).
However, as it is seen, directly deriving the voltage distribution along the axis ( ϕ i ) from this equation (i.e.Eq. ( 2)) is not feasible due to the greater number of unknowns compared to the available equations.Therefore, introducing additional conditions becomes necessary for solving the equations.
In this paper it is proposed that the axial potential can be effectively modelled by fitting it to a piecewise fifth-order spline equation.When combined with Laplace's equation, these equations can be solved together to determine the values of voltages along the axis.

Fundamental equations of the quintic spline
While a method to construct the quintic spline for a dataset with equal intervals between data points is presented in 26 , it is important to note that, in general, a lens system may have unequal gaps between the electrodes and different thicknesses for the lenses.Consequently, the set of data points created on the electrodes may have unequal intervals.As there are no references found that derive the aforementioned equations for general unequal intervals, this section addresses that gap.Therefore, in the following section, the fundamental equations of the quintic spline are obtained and formulated for the scenario in which intervals between the data points are not equal.
To initiate the formulation, let's consider a dataset represented as (x i , ϕ i ), where i = 1, 2, . . ., N , comprising N data points, each with two components.The first component is denoted as x i , and its corresponding value is ϕ i .Our goal is to fit splines to this set of N data points.Specifically, a spline is fitted between data points (x i , ϕ i ) and (x i+1 , ϕ i+1 ).
For a data set {(x i , ϕ i )}, i = 1, 2, . . ., N , the quintic spline function S(x) is defined by (3): in which S i (x) is a fifth order polynomial defined on the interval [x i , x i+1 ] and S i (x i ) = ϕ i .
To ensure a smooth transition from one spline to the next, additional conditions are imposed: where S (r) i (x) is the r th derivative of S i (x).When all these splines are combined, a generalized spline named S(x) is formed.At the data points on this generalized spline, certain conditions are imposed during spline construction, as mentioned above.
Since the polynomial is fifth order, the fourth derivative of the polynomial is a linear polynomial.Therefore the 4th derivative of S i (x) can be written as (5).
In (6), A i , B i , C i , and D i , i = 1, 2, . . ., N − 1 are coefficients which can be determined in terms of ϕ i , µ i , η i , and Z i , which are defined in (7).
Taking the derivative of (6) three times results in (8) and using the identities given by (7), results in (9), after simplification.
On substituting (21) in (18) and (24), Eqs.(25) and (26) are derived, which are referred to in this paper as the fundamental equations of the quintic spline.(2) Once the quintic spline is formulated, the axial potential has to be numerically written up to the fourth order, to be simultaneously solved together with the spline equations.
A typical rotationally symmetrical electrostatic lens system in 2D is shown in Fig. 2. For a clearer visualization of the lens geometry in 2D and 3D, along with detailed explanations of its parameters, please refer to Fig. 1.The geometry of the electrodes, including the gaps between them and the voltages at each electrode, determine the lens properties.
As aforementioned, by omitting the terms above the fourth derivatives and using the identities in (27), Eq. ( 1) can be written numerically as (28).
In electrostatic lens systems, it can be assumed that the axial potential before the first electrode and after the last electrode approaches the potential of the first and the final electrode, respectively.However, the exact position for which the above condition is valid must be identified.Referring to Fig. 2, this means the distances f 1 and f N−1 must be properly calculated.Equations ( 25) and ( 26) and ( 28) each form N − 2 equations, resulting in a total of 3N − 6 equations, while there are N unknown ϕ i , µ i , Z i for i = 1, 2, . . ., N which results in 3N unknowns in total.However, as stated above, it is assumed that ϕ 1 = V 2 and ϕ N = V N−1 , in which V 2 is the voltage of the first electrode (i.e.V EL1 ) and V N−1 is the voltage of the last electrode (i.e.V EL tot ), are both known values.Therefore, based on the above assumption, the two unknown parameters are already known, but instead, the distances f 1 and f N−1 are unknown, which again results in a total of 3N unknown parameters.Considering that there are only 3N − 6 equations, this means 6 boundary conditions must be selected first.Here we assume These assumptions are all reasonable, since before the first point and after the last point there are field-free regions and therefore assuming the second, third and fourth derivatives of voltage to be zero is also justifiable.Note that all voltages V i belonging to the same electrode j are equal to V ELj .
It should also be noted that, In Eq. ( 28), we discretized the Z axis, and as a result, the Z coordinate is embedded in the indices of the parameters.For instance, z = 0 corresponds to ϕ 1 , and with a discretization step of, for example, 0.001 mm, the point z = 0.001 is associated with ϕ 2 , and so on.When we refer to the 'first electrode, ' it implies that there are no other electrodes before it.Consequently, the voltage preceding the first electrode will be dominated by the voltage of the first electrode.This holds true for the last electrode as well.If any other part (25) of the microscope preceding the first electrode has a voltage different from that of the first electrode, we assume that part to be the first electrode.Preceding this assumed first electrode, the voltage will tend to become equal to that of the first electrode.

Calculation of f 1 and f N−1
In this section, first the values of f 1 and f N−1 are calculated for the case when there is a hole (opening in the electrode near the axis) in the first and in the last electrodes, as shown in Fig. 3a.In this case it is assumed that the voltages at the outer radius of the lens before the first electrode and after the last electrode are set to the voltage of the first electrode and the last electrode, respectively.It means that on the axis, starting from the entrance of the lens toward the left (with reference to Fig. 3a), the voltage gradually approaches the voltage of the first electrode, and at point 1 the voltage is exactly the same as that of the first electrode.Based on these conditions it can be concluded that the electric field at the first and the last point is equal to zero.The Taylor expansion of the axial potential for ϕ 2 and assuming f 1 = f can be written as: Considering the boundary conditions, and also the fact that the electric field at point 1 is zero, Eq. ( 29) can be simplified to (30) by ignoring terms with higher than fifth order derivatives.Alternatively, Eq. ( 31) can also be written.Combining (30) and (31) results in (32).
Considering Eq. ( 28) for i = 2 , Eq. ( 33) can be obtained.Equation (21) for i = 1 , by considering that in this case 1 = f and using the boundary conditions at point 1, can be written as Eq. ( 34): By combining Eqs.(33) and (34) and defining X = , Eq. ( 35) is derived: By combining (32) and (35), Eq. ( 36) is obtained.Equating and results in Eq. (37), after some simplification .( 29) (5)   1 �f 5  5! (31) ϕ (5)  1 = ϕ (4) The illustration of the lens system depicts voltages at the outer radius of the lens, with the voltages before the first electrode and after the last electrode set to the voltage of the first electrode and the last electrode, respectively.In (a), there are holes in both the first and the last electrodes.In (b), there is a hole in the first electrode, while the last electrode is closed.
Equation (37) has four roots, two of which are negative and two positive (Eq.38).The negative ones cannot be acceptable.From the two positive roots, only the smaller one (indicated in Fig. 4a) is acceptable.This is because if the larger one is selected as f 1 (as shown in Fig. 4b), considering the larger root as point 1, the smaller root falls between point 2 and point 1.In this way, there is another point on the axis that has a potential equal to V 2 .This means that the voltage before point 2 is not smoothly approaching V 2 , but oscillating, as shown in Fig. 4b, which is not possible.Therefore, the larger root cannot be an acceptable solution of the equation.
The same concept is valid for the last electrode and f N−1 can be obtained.In this case, to obtain (34), Eqs. ( 18) and ( 21) must be evaluated for i = N − 1 , and after combining them (34) is obtained.The rest of the procedure is the same and the smaller root shown in (38) can be used to calculate f N−1 with replacement of r 2 with r N−1 .
It is worth mentioning that, on the boundary conditions, if the first electrode is open (contains a hole) and the last electrode is completely closed (no hole), as shown in Fig. 3b, the final data point is ϕ N , and its voltage is known and is equal to V N , the voltage of the last electrode.As before, there are 3N − 6 equations, while there are 3N unknowns in total.However, ϕ 1 = V 2 and ϕ N = V N , which means 2 unknowns are already known, but instead, the distance f 1 must be identified, which can be found using (38).For this case it is not necessary to add another point after the last electrode and, hence, f does not exist after the last electrode.Therefore, only 5 extra boundary conditions are needed.The extra boundary conditions imposed are Note that for this case it is not necessary to set η N = 0 .In this case, while the electric field at the first point is zero, the electric field at point N is an unknown value which depends on the lens system's geometry and voltages.
Note: The term 'hole, ' as mentioned above, refers to the opening in the electrode near the axis.All electrodes must include an opening around the axis to facilitate the passage and flow of electrons, as depicted in Fig. 3a.However, the last electrode corresponds to the location of the image plane (sample plane), where it is generally closed (as shown in Fig. 3a).In some situations, it can be closed without the necessity for an opening (as shown in Fig. 3b).While such situations are infrequent, they may be of interest to electron lens designers for specific (36) www.nature.com/scientificreports/applications and therefore has been studied here as well.For a detailed exploration of these situations and their applications, refer to 12,28-30 .

Matrix formulation of the equations
First, the matrix formulations of the equations are explained for the scenario in which both the first and the last lenses have a hole in them.Assume [ϕ] N×1 is the vector of the axial potential, [µ] N×1 is the vector of the second derivative of the axial potential, and [Z] N×1 is the vector of the fourth derivative of the axial potential.Equation ( 25) in the matrix form can be written as (39).Elements of matrices Equation (26) in the matrix form can be written as (41).Elements of the matrices [C] , [D] and [K] can be obtained from (26).With the boundary conditions Z 1 = Z N = µ 1 = µ N = 0 , it is possible to reduce (41) to Eq. ( 42).In Eq. ( 42) matrices which element 1 is ϕ 1 , which is equal to the voltage of the first electrode, element N is ϕ N , which is equal to the voltage of the last electrode, and all the other elements are zero.Equation ( 28) can be written in matrix form, as in (43).In Eq. ( 43), vector [V ] represents the voltages of all points on the electrodes (therefore it is a known vector) and matrix [r] is a diagonal matrix of size (N − 2) × (N − 2) and the diagonal elements are the distances between each point on the electrode to its mapped counter point on the axis of symmetry.Equation (40) can be written as (44).Using (44) and (42), Eq. ( 45) can be calculated.Combining Eq. ( 45) and (43) results in Eq. ( 46), after simplification.Rewriting (46), Eq. (47) can be obtained from which [ϕ * ] can be calculated.Using [ϕ * ] in (45) [Z * ] is derived and using [Z * ] in (44), [µ * ] can be obtained.Adding the first and last element to the vectors [ϕ * ] , [µ * ] and [Z * ] using the boundary conditions, Z 1 = Z N = µ 1 = µ N = 0 and the fact that ϕ 1 is equal to the voltage of the first electrode and ϕ N is equal to the voltage of the last electrode, the full vectors of [ϕ] , [µ] and [Z] can be determined.

Having established [ϕ] ,
[µ] and [Z] , it is possible to formulate quintic spline equations S i (x) for i = 1, 2, • • • , N − 1 .S i (x) can be calculated using (6).Coefficients A i , B i , and D i can be directly calculated using ( 9), ( 11) and (15).To calculate C i , η i for i = 1, 2, • • • , N − 1 must first be determined from (21), and based on the assumed boundary condition η N = 0 .Having all η i and using (13), it is now possible to calculate C i .If the last electrode is closed, the same strategy as explained above can be used to calculate S i (x) The only difference is that, after determining η i for i = 1, 2, • • • , N − 1 , the value of η N should be calculated using (17) by setting i = N − 1. Note: To solve the matrix equations, various software and programming languages like MATLAB or Python can be employed for this task.In general, a linear matrix equation in the form of Ax = B can be easily solved, for example, in MATLAB by defining the matrices A and B and using the command x = A −1 B .Since all matrices presented in this Section can be computed from the geometry data and electrode voltages, the final matrices can be effortlessly calculated through a straightforward implementation in MATLAB.Following this, as explained, it is feasible to compute cubic splines, and from there, the voltage distribution along the axis of symmetry can be derived.
It is noted that for comparison with the here presented FOEM method, accurate field calculations were done using COMSOL Multiphysics, a FEM-based software tool.For mathematical formulations of solving the Poisson equation using FDM or FEM, the interested reader is referred to the literature 31-33 .

Accuracy of FOEM in field and optical parameters calculation
To analyse the accuracy of calculation of the electric field and optical parameters using FOEM, six different typical electrostatic lens systems including 3, 4 and 5 electrodes, are taken as test systems (represented in 2D in Fig. 5).The axial potential calculation is performed by the two previously existing methods of SOEM and FEM (using COMSOL 13 ) and by the newly-presented FOEM method.For the SOEM and FOEM calculation, we use MATLAB 34 .First, the axial potential and its derivatives are graphically compared.In the graphical comparison, only the graphs related to one system (system1 from Fig. 5) are depicted and presented here since all lens systems show almost similar trends.The overlaid graphs are presented in Figs. 6 and 7.
System1 : Three-electrode lens system System2: Three-electrode lens system System3: Four-electrode lens system System4: Four-electrostaƟc lens system System5: Five-electrode lens system System6: Five-electrostaƟc lens system The primary beam passing through the lens system (in red, as taken from our raytracing codes in MATLAB using FEM field calculation) is sketched for a better visualisation of the lens system.The units along the axes of the graphs are in millimetres.
In electron optical lens systems, the aberration coefficients are the factors which determine the quality of the lens systems.The lower the values of these parameters, the less the aberrations exist and therefore the better the lens system 1-3 .In our case-study, the lens systems are only suffering from the spherical and chromatic aberrations.Hence, the chromatic and spherical aberration coefficients (denoted by C c and C s ) are the deterministic factors in the lens design and optimization here.These parameters are functions of the axial potential and its derivatives.The presentation of the formulation is given in Eqs. ( 48) and (49).More details on the formulations and their derivations are skipped in this paper and can be found in the references 7,12 .
(48)   Table 1.The data related to the optical parameters such as chromatic and spherical aberrations (C c and C s ) and the image position (X c ), calculated by FEM, FOEM and SOEM, for six different typical electrostatic lens systems (presented in Fig. 5).where φ 0 is the potential at the object side, and r α (z) , refers to the principle imaging ray 1 , travelling from the object side along the optical axis, with angle 1.The principle imaging ray r α (z) is calculated by ray tracing.

Systems
The accuracy in the calculation of C S and C c is important in electron lens design.Since these parameters are functions of the potential and its derivatives, hence they can be used as an important qualitative and quantitative measure to check the accuracy of the axial potential calculations.Therefore, these optical parameters are also derived from the calculated axial potential and its derivatives by FOEM, SOEM and FEM and compared.The related graphs are presented in Fig. 8.The quantitative data related to these optical parameters are presented and compared in Tables 1 and 2.
Figure 6a, compares graphically the axial potential calculated by the three different methods (SOEM, FOEM and FEM), and a reasonably good overlap in the axial potential graphs can be seen.However, when zooming in (Fig. 6b,c, slight deviations are recognized.These deviations become more recognizable in the graph of φ ′ (Fig. 7a), and φ ′′ (Fig. 7b).It is evident that the deviations are mainly between SOEM and FEM, while there is reasonably good overlap between FOEM and FEM.
In the C s and C c graphs, shown in Fig. 8, similar to the graphs of the axial potential derivatives, the graphs of FOEM are in much higher agreement with FEM than the ones of SOEM.This behaviour is more recognizable in the C c graphs.This is because in our case-study, C s is dependent on both the first and second derivative of the axial potential, while C c is only a function of the first derivative, and the fluctuations in the second derivatives (as seen in Fig. 7a,b) are larger.
In addition to the graphical comparisons, to also achieve quantitative data comparisons, the values of optical parameters such as chromatic and spherical aberrations ( C c and C s ) as well as the image position ( X c ) (another important factor in lens design) 1 , are calculated for the six lens systems (from Fig. 5) and presented in Tables 1  and 2. The presented data of C s and C c are related to the image side.
From Table 1, it can be seen at a glance that all optical parameters, when calculated by FOEM are much closer to the accurate values calculated by FEM, compared to those calculated by SOEM.The difference between FOEM and SOEM in being more compatible with FEM, can be especially realized in C s values.
To estimate the deviation of the calculated data with FOEM and SOEM for each system, with respect to the accurate data, the error is determined by taking COMSOL data (FEM) as the reference for accurate values.The errors are presented in Table 2 for each system.The last column in Table1 presents the average value of the error for 6 lens systems.This average is computed as a mathematical mean of the errors over the six different lens systems ( ( 6 i=1 (error i ))/6)).As it is seen at a glance from Table 2, the reported errors from calculation with FOEM for all systems are much smaller than the ones calculated with SOEM.It should be noted that it might be possible to occasionally get a result with a slightly larger calculation error with FOEM than with SOEM.This scenario can occur due to the determination of the total values of C s and C c through the integration of summations involving terms dependent on φ , φ ′ and φ ′′ .Fluctuations are evident in the graphs of these parameters, particularly seen in Fig. 8a and b for phi' and phi".When examining the overlap of graphs from FOEM and SOEM with those from FEM in these figures, larger deviations are noticeable in SOEM compared to FOEM.However, in the calculations of C s and C c , owing to the complexities of the functions, the terms in the integral may occasionally compensate for each other.This compensation results in the total values of C s and C c calculated for SOEM occasionally getting closer to those calculated with FEM compared to FOEM with FEM, even though the absolute values at each point in the graphs of φ , φ ′ and φ ′′ for FOEM and FEM exhibit closer values.Such instances, however, are infrequent.In general, as shown in Table 2, the calculated data by FOEM consistently demonstrates much smaller errors compared to SOEM.
As illustrated in Table 2, FOEM values provide reasonably accurate results for X c (with the average value of error 0.7%), with respect to the accurate values of FEM.Also, noticeable is that, FOEM resulted in much smaller error percentages (approx.4 times smaller) than SOEM with respect to the accurate values.A similar trend of having higher accuracy for FOEM than SOEM is seen in the values of C c .The error percentage, on average, for C c calculation by FOEM is around 5 times smaller than by SOEM (i.e.2.8% for FOEM versus 12.79%).Having much smaller error percentage (~ 4-5 times) for C s calculation in FOEM compared to SOEM (5.6% for FOEM versus 61.6%) keeps the same trend of higher accuracy of FOEM.However, the absolute values of the error percentages in C s are much higher than in the cases of C c and X c calculation.which is due to their different functionality to the potential and its derivatives.The computation time needed to evaluate the optical parameters for each lens system, using FEM, FOEM and SOEM are ̴ 60 s, 0.5 s, and 0.4 s, respectively.
In summary, the quantitative data comparison, similar to the graphical analysis, demonstrates that FOEM performs the optical parameters calculation in a very short time with reasonably high accuracy, showing deviations ranging from around 1% to a maximum of 6% compared to FEM.Notably, FOEM exhibits significantly higher accuracy in both field calculation and optical parameter calculations compared to SOEM.This is particularly evident for optical parameters dependent on the second derivative of the potential, such as C s .

Optimization of electrostatic electron lens systems based on FOEM and a comparison with SOEM-and FEM-based optimization
As discussed earlier, the primary goal of developing a fast field calculation method for electrostatic lens systems was to incorporate it into an automated optimization routine for lens system optimization.With the identification of the proposed fast field calculation method (FOEM), the next step is to implement it in an optimization routine (here Genetic Algorithm) to assess its efficiency in producing optimized lens systems with accurately calculated (49) For the optimization, a GA is used with the generation and population of 100 and 50, respectively.Therefore, a total of 5000 systems are being evaluated in this optimization.From the results of this optimization, the 30 best non-identical systems (the system with the lowest objective function values with different designs) are taken for the comparison analysis.The related data are presented in Fig. 9.
Note: Our objective in this section is not to run a GA to assess its performance or to optimize and evaluate its effectiveness.Instead, we intend to generate several lens systems, surpassing the six systems listed in Table 1, each with distinct geometries.The primary goal is to reevaluate the accuracy of calculations for C S , C C , and X C by FOEM compared to FEM and SOEM across a broader range of lens systems.Additionally, we aim to investigate how the trend of decreasing C S and C C , observed in FOEM calculations, translates when the same lens systems are computed using FEM.Specifically, we want to understand whether C S and C C maintain a similar trend of decrement when they are calculated by FOEM (compared to when they are calculated by FEM).This evaluation serves two purposes: firstly, to confirm the accuracy of C S , C C , and X C calculations for various systems, not just those in Table 1, and secondly, to find out whether lens systems with identical X C values with different C S and C C values when calculated by FOEM, would show a similar trend of decrement of optical parameters when using FEM.If that can be established, it suggests that FOEM is a potentially valuable method to be implemented in optimization routines for optical parameter calculations, and it would guide us towards systems with a decrement in their C S and C C , even if there are some deviations in the accurately calculated values of C S and C C .Simultaneously, a similar comparison is also performed for SOEM, allowing us to compare the results when they were calculated by FOEM compared to SOEM.To conduct this analysis, we could have manually found 30 diverse lens systems, each with X C at 15 mm.However, the task to manually create these test lens systems is challenging.To address this, instead, we utilized a GA run, incorporating a constraint of X C = 15 mm .Subsequently, we selected 30 lens systems that met the constraint and presented varying geometries.This automated approach allowed us to efficiently find 30 systems with decreasing C S and C C values, facilitating our evaluation.We didn't include all 5000 systems from their respective generations in the graph since they were too crowded.Instead, we manually selected the 30 best non-identical systems-those with the lowest objective function values and different designs-for this comparison.The reason for not selecting the real 30 best systems out of the GA run was to avoid having almost 30 identical systems that GA converged to as solutions in the last generations.This would have hindered a proper observation of the trend of decrement.Additionally, not having 30 almost very different lens systems would not provide a suitable basis for testing optical parameter accuracy calculations.www.nature.com/scientificreports/Therefore, we manually chose 30 systems to represent the overall trend of decrement in their objective functions.These systems were sorted in the plot and identified by their sorted system numbers (System 1-30) in Fig. 8.The data of X C , C C and C S , for these 30 systems are given as squares (in red) in Fig. 9a-c, respectively.The optical parameters of X C , C C and C S for these 30 systems are again calculated using SOEM and FEM (by COMSOL).The data calculated by SOEM are presented with stars (in blue), and by FEM are shown with circles (in black) in the same Fig.9a-d).As it is clear from Fig. 9a, the data for X C are reasonably accurate both for FOEM and SOEM compared to FEM, but with higher accuracy for FOEM than for SOEM.It is noticeable from Fig. 9b-d, that the data of C C and C S for FEM and SOEM follow the same decreasing trend as FOEM, however, with much less deviation from the accurate data of FEM for FOEM data compared to those of SOEM.
The first conclusion is that both SOEM and FOEM can be used in an initial optimization to start from no initial data and reach reasonably optimized systems.However, it is important to note that the inaccuracy in the data obtained with SOEM is large, specifically in C S (see Fig. 9c and zoomed-in in Fig. 9d).

The efficiency of GA optimization based on FOEM, FEM and SOEM
In this subsection, the objective is to compare the outcomes of FOEM + GA, FEM + GA, and SOEM + GA.To achieve this, different optimizations, as detailed below, will be executed and subsequently compared: FEM + GA optimization starting with random parameters (very time consuming), then FEM + GA starting with a feed from SOEM and finally FEM + GA starting with a feed from FOEM.The generation (denoted by 'G') and population (denoted by 'P') in GA are taken as 100 and 50, respectively.The results of these three optimizations are presented in Fig. 10a, 11a and 12a.Figures 10b, 11b and 12b depict a zoomed-in view.As previously discussed, the results from FEM are accurate.Therefore, the results taken from "FEM + GA" do not need to be fine-tuned.In order to get accurate and optimized out-put results, the 10 best optimized systems are taken and fed into another optimization by GA as its initial population, but this time the objective function is calculated using FEM.The results of these optimizations which are called "SOEM + FEM + GA" and "FOEM + FEM + GA", are presented in Figs.11c and 12c, respectively.The quantitative data related to these optimizations are represented and summarised in Table 3.
As is seen from Fig. 10a, to run the full optimization for "FEM + GA", the computation time (denoted by CT in the Fig. 10, 11 and 12), is 80 h.However, the best system, with the objective function of 15.75 (shown by ' A' in Fig. 10b), is reached at generation G = 83.Hence, the computation time needed to reach the best system in this case is around 70 h.After this time, the optimization could not be improved any further within its 100 generations.Therefore the CT needed for "FEM + GA" can be reported as 70 h and the best system has OF = 15.75 (represented in Table 3 as well).
Figure 11, shows the results of "SOEM + FEM + GA".The first plot, 9a, presents the improvement of the objective function versus generation during the optimization.The total computation time to run this optimization is CT ~ 33 min.The zoomed-in data is shown in Fig. 11b where it is seen that the optimization has improvements almost until the last generation.The best system (system B in Fig. 11b) is found at generation of G = 99, with an objective function of OF = 19.01.The computation time to find this system is ~ 32 min.The related data for this system is given in Table 3.As is evident from Table 3, the image position of this system when derived from "SOEM + GA" is reported to be X C _SOEM = 14.96 mm.However, when the optical parameters of this system are calculated accurately by FEM, X C _FEM = 14.54 mm and OF_FEM = 13.56 is achieved.optimized system (system C' in Fig. 12c) after only the first generation (at G = 2, OF = 14.32).After the second generation the optimization improvements are almost negligible.This means that the initial systems which had been fed were close to the accurate solutions calculated by FEM.
The reason why the "SOEM + FEM + GA" optimization takes a longer time than "FOEM + FEM + GA" to find the optimized system is that the objective function and image position of the optimized system found by this optimization (system B) deviate more from the corresponding accurately calculated values by FEM (see Table 3) than those achieved by "FOEM + GA" (system C).
In summary, based on what has been discussed above and reported in Table 3, it is seen that the total computation time needed to run a fully-automated optimization, for "FEM + GA", "SOEM + FEM + GA" and "FOEM + FEM + GA", is 70 h, 7 h and 1.25 h (i.e.35 min + 40 min), respectively.The optimized systems found by these three methods, which satisfied the constraints, have objective function values of 15.75, 15.66 and 14.32, respectively.Although "SOEM + FEM + GA" could reach an optimized system with low objective function within a reasonably short time (~ 7 h), "FOEM + FEM + GA" could reach it ̴ 10 times faster and it is even ̴ 70 times faster than the optimization based on "FEM + GA".From these results, it is clear that "FOEM + FEM + GA" could reach a lower value of the objective function and in an enormously shorter computation time than the other two optimization strategies.
However, what is noticeable here is that, as it is seen from the last row of the Table 3, the optimized system found by only "FOEM + GA" got close enough to the accurate optimized system, which after a very short time of tuning (10 min) could achieve an optimized system which could not be further improved by a full optimization .System B is the best system found in this optimization.'CT' stands for the computation time needed to run the optimization until the generation indicated by the green arrow.Figure (c) presents the result of an optimization using GA, called "FOEM + FEM + GA", with the initial data fed from a previous optimization of "FOEM + GA".A lens system depicted by C' , is an optimized system resulting from "FOEM + FEM + GA" with OF = 14.32, occurring at generation 2.
Hence, it is evident that FOEM results are close enough to the accurate results such that this method has the possibility to be even used separately in an optimization run without the need of another full run of FEM, while SOEM needed that.This method therefore can have high impact on the process of automated design and optimization of electrostatic lens systems.
To have a visualisation of the electrostatic lens system optimized by "FOEM + GA", system C' is presented in Fig. 13 (2D in a, 3D in b).

Conclusion
In this work, a novel fast field calculation method, with reasonably high accuracy, for cylindrical rotationallysymmetrical electrostatic lens systems is presented.It is performed by solving the Laplace equation while keeping the terms up to its fourth order derivatives and utilizing quintic splines.This method is therefore called Fourth Order Electrode Method (FOEM) by the authors.First, fundamental equations of the quintic spline are derived for a system with unequal distances between the points, which is new as far as we know.Using the Laplace equation and by means of discretizing the axis-of-symmetry and approximation of the axial potential in each Table 3.The data related to the comparison between three optimization strategies of "FEM + GA", "FOEM + FEM + GA" and "SOEM + FEM + GA".

Figure 1 .
Figure 1.(a) Illustration of the rotationally symmetrical electrostatic lens system with 4 electrodes in a 3D representation.(b) Schematic depicting the rotationally symmetrical electrostatic lens system with 4 electrodes in a 2D representation.Here, T j , R j , and V EL j denote the thicknesses, radii, and voltages at each electrode.G j represents the gaps between two consecutive electrodes.

Figure 2 .
Figure 2. A typical rotationally symmetrical electrostatic lens system in 2D.The z-axis is the axis of rotational symmetry.

Figure 4 .
Figure 4. Proper selection of f 1 and f N−1 from the roots of (38).(a) axial potential smoothly approaches the voltage of the first electrode when f is equal to the smaller root of (38), (b) axial potential oscillates before it reaches the voltage of the first electrode when f is equal to the larger root of (38).
[A] and [B] can be obtained from (25).The close format of matrices are represented in the text.The open-format of matrices are given in the appendix for more clarity.With the boundary conditions Z 1 = Z N = µ 1 = µ N = 0 , it is possible to reduce Eq.(39) and obtain (40).In Eq. (40) [Z * ] and [µ * ] are the vectors of the fourth derivative and second derivative of axial potential, respectively, starting from point 2 until point N − 1 .Matrices [A * ] and [B * ] are the reduced [A] and [B].

Figure 5 .
Figure 5. Six different typical electrostatic lens systems (in 2D) taken as test systems for making the comparison of the methods of SOEM, FOEM and FEM.The colours indicate the voltages at each electrode.The primary beam passing through the lens system (in red, as taken from our raytracing codes in MATLAB using FEM field calculation) is sketched for a better visualisation of the lens system.The units along the axes of the graphs are in millimetres.

Figure 6 .
Figure 6.(a) The overlapped graphs of the axial potential for System 1 from Fig. 5, calculated with different methods of SOEM, FOEM and FEM.(b) and (c) the enlarged sections of the graphs.

Figure 7 .
Figure 7.The overlapped graphs of the first (a) and second (b) derivative of the axial potential for System1 from Fig. 5, calculated with different methods of SOEM, FOEM and FEM.

Figure 8 .
Figure 8.The overlapped graphs of chromatic (a) and spherical (b) aberration coefficients for System1 from Fig. 5, calculated with different methods of SOEM, FOEM and FEM.The unit of the y-axis is in millimeters.

Figure 9 .
Figure 9. (a) Comparison of crossover point ( X C ), (b) chromatic aberation cofficient ( C C ), (c,d) and spherical aberration coefficient ( C S ), for the best 30 systems obtained using SOEM, FOEM and COMSOL.The 30 systems are the best non-identical optimized systems derived from an optimization using GA with data calculation by FOEM.

Figure 10 .
Figure 10.(a) Optimization of a 4-electrode lens system using GA, while the optical parameters are calculated by FEM ("FEM + GA").The graph represents the objective function (OF) values (Y-axis) versus generation (X-axis).The black/blue dots show the best/mean objective function values at each generation.Figure (b) shows the zoomed-in data of figure a. System A is the best system found in this optimization with OF = 15.75, occurring at generation 83.'CT' stands for the computation time needed to run the optimization till the generation indicated by the green arrow.

Figure 12 .
Figure 12.(a) The results from an optimization of a 4-electrostatic lens system using GA, where the optical parameters are calculated by FOEM ("FOEM + GA").The graph represents the objective function (OF) values (Y-axis) versus generation (X-axis).The black/blue dots show the best/mean objective function values at each generation.Figure (b) shows the zoomed-in data of figure (a).System B is the best system found in this optimization.'CT' stands for the computation time needed to run the optimization until the generation indicated by the green arrow.Figure(c) presents the result of an optimization using GA, called "FOEM + FEM + GA", with the initial data fed from a previous optimization of "FOEM + GA".A lens system depicted by C' , is an optimized system resulting from "FOEM + FEM + GA" with OF = 14.32, occurring at generation 2.

Table 2 .
The percentage of error in calculation of the optical parameter derived by FOEM and SOEM, for six different typical electrostatic lens systems (presented in Fig.5) compared with the data calculated with the accurate method of FEM.