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, \ldots ,tot\), where \(tot\) is the total number of electrodes). While \(G_{j}\) represents the gaps between two consecutive electrodes ((\(j = 1,2, \ldots ,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}

$$\varphi \left( {r,z} \right) = \varphi \left( {0,z} \right) – \frac{{r^{2} }}{4}\varphi^{\left( 2 \right)} \left( {0,z} \right) + \frac{{r^{4} }}{64}\varphi^{\left( 4 \right)} \left( {0,z} \right) + \cdots + \frac{{\left( { – 1} \right)^{m} r^{2m} }}{{\left( {m!} \right)^{2} 2^{2m} }}\varphi^{{\left( {2m} \right)}} \left( {0,z} \right) + \cdots$$

(1)

in which \(r\) is the radial distance to the axis, \(z\) is the horizontal distance from the entrance of the lens system, \(\varphi \left( {r,z} \right)\) is the potential in space at coordinate \(\left( {r,z} \right)\), \(\varphi \left( {0,z} \right)\) is the axial potential at coordinate \(\left( {0,z} \right)\), and \(\varphi^{{\left( {2i} \right)}} \left( {0,z} \right)\) is the \(\left( {2i} \right)\)th derivative of the axial potential at coordinate \(\left( {0,z} \right)\). 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:

$$V_{i} \left( {r_{i} } \right) = \varphi_{i} – \frac{{r_{i}^{2} }}{4}\varphi_{i}^{\left( 2 \right)} + \frac{{r_{i}^{4} }}{64}\varphi_{i}^{\left( 4 \right)} \quad i = 2, 3, \ldots , N – 1$$

(2)

Considering Eq.Â (2), it is evident that the potential distribution along the axis (\(\varphi_{i} )\) is linked to voltages at the surface of the electrodes (\(V_{i} \left( {r_{i} } \right))\), where the values of \(V_{i} \left( {r_{i} } \right))\) are known at each electrode and correspond to the electrode voltages assigned to that specific electrode (\(V_{ELj} { }, j = 1,2, \ldots ,tot\)).

However, as it is seen, directly deriving the voltage distribution along the axis (\(\varphi_{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 \(\left( {x_{i} ,\varphi_{i} } \right),\) where \({ }i = 1,2, \ldots ,{ }N\), comprising \(N\) data points, each with two components. The first component is denoted as \(x_{i}\), and its corresponding value is \(\varphi_{i}\). Our goal is to fit splines to this set of \(N\) data points. Specifically, a spline is fitted between data points \(\left( {x_{i} ,\varphi_{i} } \right){ }\) and \(\left( {x_{i + 1} ,\varphi_{i + 1} } \right).\)

For a data set \(\left\{ {\left( {x_{i} ,\varphi_{i} } \right)} \right\}, i = 1, 2, \ldots , N\), the quintic spline function \(S\left( x \right)\) is defined by (3):

$$S\left( x \right) = S_{i} \left( x \right), \quad x \in \left[ {x_{i} , x_{i + 1} } \right], \;\;i = 1, 2, \ldots , N – 1$$

(3)

in which \(S_{i} \left( x \right)\) is a fifth order polynomial defined on the interval \(\left[ {x_{i} , x_{i + 1} } \right]\) and \(S_{i} \left( {x_{i} } \right) = \varphi_{i}\).

To ensure a smooth transition from one spline to the next, additional conditions are imposed:

$$S_{i}^{\left( r \right)} \left( {x_{i + 1} } \right) = S_{i + 1}^{\left( r \right)} \left( {x_{i + 1} } \right),\quad r = 0,{ }1,{ }2,{ }3,{ }4{ }$$

(4)

where \(S_{i}^{\left( r \right)} \left( x \right)\) is the \(r\)th derivative of \(S_{i} \left( x \right)\).

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} \left( x \right)\) can be written as (5).

$$S_{i}^{\left( 4 \right)} \left( x \right) = Z_{i + 1} \frac{{\left( {x – x_{i} } \right)}}{{\Delta_{i} }} + Z_{i} \frac{{\left( {x_{i + 1} – x} \right)}}{{\Delta_{i} }}$$

(5)

in which \(Z_{i} = S_{i}^{\left( 4 \right)} \left( {x_{i} } \right),\;\; x \in \left[ {x_{i} , x_{i + 1} } \right]\), and \(\Delta_{i} = x_{i + 1} – x_{i}\). Integrating (5) four times, and after some rearrangements, Eq.Â (6) can be obtained.

$$S_{i} \left( x \right) = Z_{i + 1} \frac{{\left( {x – x_{i} } \right)^{5} }}{{120\Delta_{i} }} + Z_{i} \frac{{\left( {x_{i + 1} – x} \right)^{5} }}{{120\Delta_{i} }} + A_{i} \left( {x – x_{i} } \right)^{3} + B_{i} \left( {x_{i + 1} – x} \right)^{2} + C_{i} \left( {x – x_{i} } \right) + D_{i} \left( {x_{i + 1} – x} \right)$$

(6)

In (6), \(A_{i} , B_{i} , C_{i} ,\) and \(D_{i}\), \(i = 1, 2, \ldots , N – 1\) are coefficients which can be determined in terms of \(\varphi_{i} , \mu_{i} , \eta_{i} ,\) and \(Z_{i}\), which are defined in (7).

$$\varphi_{i} = S_{i} \left( {x_{i} } \right) ,\;\; \mu_{i} = S_{i}^{\left( 2 \right)} \left( {x_{i} } \right) , \;\;\eta_{i} = S_{i}^{\left( 3 \right)} \left( {x_{i} } \right) ,\;\;Z_{i} = S_{i}^{\left( 4 \right)} \left( {x_{i} } \right)$$

(7)

Taking the derivative of (6) three times results in (8) and using the identities given by (7), results in (9), after simplification.

$$S_{i}^{\left( 3 \right)} \left( x \right) = Z_{i + 1} \frac{{\left( {x – x_{i} } \right)^{2} }}{{2\Delta_{i} }} – Z_{i} \frac{{\left( {x_{i + 1} – x} \right)^{2} }}{{2\Delta_{i} }} + 6A_{i}$$

(8)

$$A_{i} = \frac{{\Delta_{i} }}{12}Z_{i} + \frac{{\eta_{i} }}{6}$$

(9)

EquationÂ (10) is then obtained by taking the derivative of Eq.Â (6) twice, which produces Eq.Â (11), using the identities given by Eq.Â (7) and upon simplification.

$$S_{i}^{\left( 2 \right)} \left( x \right) = Z_{i + 1} \frac{{\left( {x – x_{i} } \right)^{3} }}{{6\Delta_{i} }} + Z_{i} \frac{{\left( {x_{i + 1} – x} \right)^{3} }}{{6\Delta_{i} }} + 6A_{i} \left( {x – x_{i} } \right) + 2B_{i}$$

(10)

$$B_{i} = \frac{{\mu_{i} }}{2} – \frac{{\Delta_{i}^{2} }}{12}Z_{i}$$

(11)

Replacing \(x = x_{i + 1}\) in (6) and using the identities given by (7) yields (12) which produces (13), after simplification.

$$\varphi_{i + 1} = Z_{i + 1} \frac{{\Delta_{i}^{4} }}{120} + A_{i} \Delta_{i}^{3} + C_{i} \Delta_{i}$$

(12)

$$C_{i} = \frac{{\varphi_{i + 1} }}{{\Delta_{i} }} – \frac{{\Delta_{i}^{3} }}{12}Z_{i} – \frac{{\Delta_{i}^{3} }}{120}Z_{i + 1} – \frac{{\Delta_{i}^{2} }}{6}\eta_{i}$$

(13)

Equation (14) is obtained by replacing \(x = x_{i}\) in (6) and using the identities given by (7). Following simplification, using (11) in (14) then produces (15).

$$D_{i} = \frac{{\varphi_{i} }}{{\Delta_{i} }} – Z_{i} \frac{{\Delta_{i}^{3} }}{120} – B_{i} \Delta_{i}$$

(14)

$$D_{i} = \frac{{\varphi_{i} }}{{\Delta_{i} }} + \frac{{9\Delta_{i}^{3} }}{120}Z_{i} – \frac{{\mu_{i} }}{2}\Delta_{i}$$

(15)

In the equations below, the conditions of continuity of derivatives given in (4) for \(r = 1, 2, 3\) are applied. For the case when \(r = 3\), after simple evaluation, (16) and (17) are obtained. By equating (16) and (17) and using (9), Eq.Â (18) is derived.

$$S_{i}^{\left( 3 \right)} \left( {x_{i + 1} } \right) = Z_{i + 1} \frac{{\Delta_{i} }}{2} + 6A_{i}$$

(16)

$$S_{i + 1}^{\left( 3 \right)} \left( {x_{i + 1} } \right) = – Z_{i + 1} \frac{{\Delta_{i + 1} }}{2} + 6A_{i + 1}$$

(17)

$$\frac{{\Delta_{i} }}{2}Z_{i + 1} + \frac{{\Delta_{i} }}{2}Z_{i} + \eta_{i} – \eta_{i + 1} = 0$$

(18)

For the case when \(r = 2,\) through a straightforward assessment, (19) and (20) are obtained. By equating (19) and (20) and using (9) and (11), Eq.Â (21) is produced.

$$S_{i}^{\left( 2 \right)} \left( {x_{i + 1} } \right) = Z_{i + 1} \frac{{\Delta_{i}^{2} }}{6} + 6A_{i} \Delta_{i} + 2B_{i}$$

(19)

$$S_{i + 1}^{\left( 2 \right)} \left( {x_{i + 1} } \right) = Z_{i + 1} \frac{{\Delta_{i + 1}^{2} }}{6} + 2B_{i + 1}$$

(20)

$$\eta_{i} = – \frac{{2\Delta_{i} }}{6}Z_{i} – \frac{{\Delta_{i} }}{6}Z_{i + 1} – \frac{{\mu_{i} }}{{\Delta_{i} }} + \frac{{\mu_{i + 1} }}{{\Delta_{i} }}$$

(21)

For the case when \(r = 1\) and after simple evaluation, (22) and (23) are obtained. By equating (22) and (23) and using (9), (11), (13), and (15), Eq.Â (24) then results.

$$S_{i}^{\left( 1 \right)} \left( {x_{i + 1} } \right) = Z_{i + 1} \frac{{\Delta_{i}^{3} }}{24} + 3A_{i} \Delta_{i}^{2} + C_{i} – D_{i}$$

(22)

$$S_{i + 1}^{\left( 1 \right)} \left( x \right) = – Z_{i + 1} \frac{{\Delta_{i + 1}^{3} }}{24} – 2B_{i + 1} \Delta_{i + 1} + C_{i + 1} – D_{i + 1}$$

(23)

$$\frac{{11\Delta_{i}^{3} }}{120}Z_{i} + \frac{{4\Delta_{i}^{3} + 4\Delta_{i + 1}^{3} }}{120}Z_{i + 1} + \frac{{\Delta_{i + 1}^{3} }}{120}Z_{i + 2} + \frac{{2\Delta_{i}^{2} }}{6}\eta_{i} + \frac{{\Delta_{i + 1}^{2} }}{6}\eta_{i + 1} + \frac{{\Delta_{i} }}{2}\mu_{i} + \frac{{\Delta_{i + 1} }}{2}\mu_{i + 1} = \frac{{\varphi_{i + 2} }}{{\Delta_{i + 1} }} – \frac{{\varphi_{i + 1} }}{{\Delta_{i + 1} }} – \frac{{\varphi_{i + 1} }}{{\Delta_{i} }} + \frac{{\varphi_{i} }}{{\Delta_{i} }}$$

(24)

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.

$$\frac{{\Delta_{i} }}{6}Z_{i} + \frac{{2(\Delta_{i} + \Delta_{i + 1} )}}{6}Z_{i + 1} + \frac{{\Delta_{i + 1} }}{6}Z_{i + 2} – \frac{{\mu_{i} }}{{\Delta_{i} }} + \frac{{\mu_{i + 1} }}{{\Delta_{i} }} + \frac{{\mu_{i + 1} }}{{\Delta_{i + 1} }} – \frac{{\mu_{i + 2} }}{{\Delta_{i + 1} }} = 0\quad i = 1, 2, \ldots , N – 2$$

(25)

$$- \frac{{7\Delta_{i}^{3} }}{360}Z_{i} – \frac{{8\left( {\Delta_{i}^{3} + \Delta_{i + 1}^{3} } \right)}}{360}Z_{i + 1} – \frac{{7\Delta_{i + 1}^{3} }}{360}Z_{i + 2} + \frac{{\Delta_{i} }}{6}\mu_{i} + \frac{{2\left( {\Delta_{i} + \Delta_{i + 1} } \right)}}{6}\mu_{i + 1} + \frac{{\Delta_{i + 1} }}{6}\mu_{i + 2} = \frac{{\varphi_{i + 2} }}{{\Delta_{i + 1} }} – \frac{{\varphi_{i + 1} }}{{\Delta_{i + 1} }} – \frac{{\varphi_{i + 1} }}{{\Delta_{i} }} + \frac{{\varphi_{i} }}{{\Delta_{i} }} \quad i = 1, 2, \ldots , N – 2$$

(26)

### Solving the Laplace equations up to the fourth order derivatives

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

$$\mu_{i} = \varphi_{i}^{\left( 2 \right)} \left( {0,z} \right) ,\;\;Z_{i} = \varphi_{i}^{\left( 4 \right)} \left( {0,z} \right)$$

(27)

$$V_{i} \left( {r_{i} } \right) = \varphi_{i} – \frac{{r_{i}^{2} }}{4}\mu_{i} + \frac{{r_{i}^{4} }}{64}Z_{i} \quad i = 2, 3, \ldots , N – 1$$

(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 \(\Delta f_{1}\) and \(\Delta 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 \(\varphi_{i}\), \(\mu_{i} , Z_{i}\) for \(i = 1, 2, \ldots ,N\) which results in \(3N\) unknowns in total. However, as stated above, it is assumed that \(\varphi_{1} = V_{2}\) and \(\varphi_{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 \(\Delta f_{1}\) and \(\Delta 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 \(\mu_{1} = \mu_{N} = \eta_{1} = \eta_{N} = Z_{1} = Z_{N} = 0\). 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 \(\varphi_{1}\), and with a discretization step of, for example, 0.001Â mm, the point zâ=â0.001 is associated with \(\varphi_{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 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 \(\Delta {\varvec{f}}_{1}\) and \(\Delta {\varvec{f}}_{{{\varvec{N}} – 1}}\)

In this section, first the values of \(\Delta f_{1}\) and \(\Delta 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 \(\varphi_{2}\) and assuming \(\Delta f_{1} = \Delta f\) can be written as:

$$\varphi_{2} = \varphi_{1} + \varphi_{1}{\prime} \frac{\Delta f}{{1!}} + \varphi_{1}^{^{\prime\prime}} \frac{{\Delta f^{2} }}{2!} + \varphi_{1}^{\left( 3 \right)} \frac{{\Delta f^{3} }}{3!} + \varphi_{1}^{\left( 4 \right)} \frac{{\Delta f^{4} }}{4!} + \varphi_{1}^{\left( 5 \right)} \frac{{\Delta f^{5} }}{5!} + \cdots$$

(29)

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

$$\varphi_{2} = \varphi_{1} + \varphi_{1}^{\left( 5 \right)} \frac{{\Delta f^{5} }}{5!}$$

(30)

$$\varphi_{1}^{\left( 5 \right)} = \frac{{\varphi_{2}^{\left( 4 \right)} – \varphi_{1}^{\left( 4 \right)} }}{\Delta f} = \frac{{\varphi_{2}^{\left( 4 \right)} }}{\Delta f}$$

(31)

$$\varphi_{2} = \varphi_{1} + \varphi_{2}^{\left( 4 \right)} \frac{{\Delta f^{4} }}{5!}$$

(32)

Considering Eq.Â (28) for \(i = 2\), Eq.Â (33) can be obtained.

$$V_{2} = \varphi_{2} – \frac{{r_{2}^{2} }}{4}\varphi_{2}^{\left( 2 \right)} + \frac{{r_{2}^{4} }}{64}\varphi_{2}^{\left( 4 \right)}$$

(33)

EquationÂ (21) for \(i = 1\), by considering that in this case \(\Delta_{1} = \Delta f\) and using the boundary conditions at point 1, can be written as Eq.Â (34):

$$\varphi_{2}^{\left( 2 \right)} = \frac{{\left( {\Delta f} \right)^{2} }}{6}\varphi_{2}^{\left( 4 \right)}$$

(34)

By combining Eqs.Â (33) and (34) and defining \(X = \left( {\frac{{r_{2}^{4} }}{64} – \frac{{r_{2}^{2} }}{4}\frac{{\left( {\Delta f} \right)^{2} }}{6}} \right)\), Eq.Â (35) is derived:

$$\varphi_{2}^{\left( 4 \right)} = \frac{{V_{2} – \varphi_{2} }}{X}$$

(35)

By combining (32) and (35), Eq.Â (36) is obtained. Equating \(\frac{{\Delta f^{4} }}{5!X} = – 1\) produces the desired result of \(\varphi_{1} = V_{2}\), and results in Eq.Â (37), after some simplification .

$$\varphi_{1} = \varphi_{2} \left( {1 + \frac{{\Delta f^{4} }}{5!X}} \right) – V_{2} \frac{{\Delta f^{4} }}{5!X}$$

(36)

$$\Delta f^{4} – 5r_{2}^{2} \Delta f^{2} + \frac{15}{8}r_{2}^{4} = 0$$

(37)

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 \(\Delta 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.

$$\Delta f = r_{2} \sqrt {\frac{{5 – \sqrt{\frac{35}{2}} }}{2}} , \;\;\Delta f = r_{2} \sqrt {\frac{{5 + \sqrt{\frac{35}{2}} }}{2}}$$

(38)

The same concept is valid for the last electrode and \(\Delta 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 \(\Delta 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 \(\varphi_{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, \(\varphi_{1} = V_{2}\) and \(\varphi_{N} = V_{N}\), which means 2 unknowns are already known, but instead, the distance \(\Delta 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, \(\Delta f\) does not exist after the last electrode. Therefore, only 5 extra boundary conditions are needed. The extra boundary conditions imposed are \(\mu_{1} = \mu_{N} = \eta_{1} = Z_{1} = Z_{N} = 0\). Note that for this case it is not necessary to set \(\eta_{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 applications and therefore has been studied here as well. For a detailed exploration of these situations and their applications, refer to^{12,28,29,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 \(\left[ \varphi \right]_{N \times 1}\) is the vector of the axial potential, \(\left[ \mu \right]_{N \times 1}\) is the vector of the second derivative of the axial potential, and \(\left[ Z \right]_{N \times 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 \(\left[ A \right]\) and \(\left[ B \right]\) 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.

$$\left[ A \right]_{{\left( {N – 2} \right) \times N}} \left[ Z \right]_{N \times 1} = \left[ B \right]_{{\left( {N – 2} \right) \times N}} \left[ \mu \right]_{N \times 1}$$

(39)

With the boundary conditions \(Z_{1} = Z_{N} = \mu_{1} = \mu_{N} = 0\), it is possible to reduce Eq.Â (39) and obtain (40). In Eq.Â (40) \(\left[ {Z^{*} } \right]\) and \(\left[ {\mu^{*} } \right]\) are the vectors of the fourth derivative and second derivative of axial potential, respectively, starting from point 2 until point \(N – 1\). Matrices \(\left[ {A^{*} } \right]\) and \(\left[ {B^{*} } \right]\) are the reduced \(\left[ A \right]\) and \(\left[ B \right]\).

$$\left[ {A^{*} } \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}} \left[ {Z^{*} } \right]_{{\left( {N – 2} \right) \times 1}} = \left[ {B^{*} } \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}} \left[ {\mu^{*} } \right]_{{\left( {N – 2} \right) \times 1}}$$

(40)

EquationÂ (26) in the matrix form can be written as (41). Elements of the matrices \(\left[ C \right]\), \(\left[ D \right]\) and \(\left[ K \right]\) can be obtained from (26). With the boundary conditions \(Z_{1} = Z_{N} = \mu_{1} = \mu_{N} = 0\), it is possible to reduce (41) to Eq.Â (42). In Eq.Â (42) matrices \(\left[ {C^{*} } \right]\), \(\left[ {D^{*} } \right]\) and \(\left[ {K^{*} } \right]\) are the reduced \(\left[ C \right]\), \(\left[ D \right]\) and \(\left[ K \right]\). In (42) \(\left[ {U_{0} } \right]\) is a vector of size \(\left( {N – 2} \right) \times 1\), in which element 1 is \(\varphi_{1}\), which is equal to the voltage of the first electrode, element N is \(\varphi_{N}\), which is equal to the voltage of the last electrode, and all the other elements are zero.

$$\left[ C \right]_{{\left( {N – 2} \right) \times N}} \left[ Z \right]_{N \times 1} + \left[ D \right]_{{\left( {N – 2} \right) \times N}} \left[ \mu \right]_{N \times 1} = \left[ K \right]_{{\left( {N – 2} \right) \times N}} \left[ \varphi \right]_{N \times 1}$$

(41)

$$\left[ {C^{*} } \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}} \left[ {Z^{*} } \right]_{{\left( {N – 2} \right) \times 1}} + \left[ {D^{*} } \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}} \left[ {\mu^{*} } \right]_{{\left( {N – 2} \right) \times 1}} = \left[ {K^{*} } \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}} \left[ {\varphi^{*} } \right]_{{\left( {N – 2} \right) \times 1}} + \left[ {U_{0} } \right]_{{\left( {N – 2} \right) \times 1}}$$

(42)

EquationÂ (28) can be written in matrix form, as in (43). In Eq.Â (43), vector \(\left[ V \right]\) represents the voltages of all points on the electrodes (therefore it is a known vector) and matrix \(\left[ r \right]\) is a diagonal matrix of size \(\left( {N – 2} \right) \times \left( {N – 2} \right)\) and the diagonal elements are the distances between each point on the electrode to its mapped counter point on the axis of symmetry.

$$\left[ V \right]_{{\left( {N – 2} \right) \times 1}} = \left[ {\varphi^{*} } \right]_{{\left( {N – 2} \right) \times 1}} – \frac{{\left[ r \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}}^{2} }}{4}\left[ {\mu^{*} } \right]_{{\left( {N – 2} \right) \times 1}} + \frac{{\left[ r \right]_{{\left( {N – 2} \right) \times \left( {N – 2} \right)}}^{4} }}{64}\left[ {Z^{*} } \right]_{{\left( {N – 2} \right) \times 1}}$$

(43)

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 \(\left[ {\varphi^{*} } \right]\) can be calculated. Using \(\left[ {\varphi^{*} } \right]\) in (45) \(\left[ {Z^{*} } \right]\) is derived and using \(\left[ {Z^{*} } \right]\) in (44), \(\left[ {\mu^{*} } \right]\) can be obtained. Adding the first and last element to the vectors \(\left[ {\varphi^{*} } \right]\), \(\left[ {\mu^{*} } \right]\) and \(\left[ {Z^{*} } \right]\) using the boundary conditions, \(Z_{1} = Z_{N} = \mu_{1} = \mu_{N} = 0\) and the fact that \(\varphi_{1}\) is equal to the voltage of the first electrode and \(\varphi_{N}\) is equal to the voltage of the last electrode, the full vectors of \(\left[ \varphi \right]\), \(\left[ \mu \right]\) and \(\left[ Z \right]\) can be determined.

$$\left[ {\mu^{*} } \right] = \left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]\left[ {Z^{*} } \right]$$

(44)

$$\left[ {Z^{*} } \right] = \left\{ {\left[ C \right] + \left[ D \right]\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right\}^{ – 1} \left( {\left[ {K^{*} } \right]\left[ {\varphi^{*} } \right] + \left[ {U_{0} } \right]} \right)$$

(45)

$$\begin{aligned} & \left[ V \right] – \left( {\frac{{\left[ r \right]^{4} }}{64} – \frac{{\left[ r \right]^{2} }}{4}\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right)\left\{ {\left[ C \right] + \left[ D \right]\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right\}^{ – 1} \left[ {U_{0} } \right] \\ & \quad = \left( {I + \left( {\frac{{\left[ r \right]^{4} }}{64} – \frac{{\left[ r \right]^{2} }}{4}\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right)\left\{ {\left[ C \right] + \left[ D \right]\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right\}^{ – 1} \left[ {K^{*} } \right]} \right)\left[ {\varphi^{*} } \right] \\ \end{aligned}$$

(46)

$$\begin{aligned} & \left[ {\varphi^{*} } \right] = \left[ P \right]^{ – 1} \left[ Q \right] \\ & \left[ P \right] = \left( {I + \left( {\frac{{\left[ r \right]^{4} }}{64} – \frac{{\left[ r \right]^{2} }}{4}\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right)\left\{ {\left[ C \right] + \left[ D \right]\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right\}^{ – 1} \left[ {K^{*} } \right]} \right) \\ & \left[ Q \right] = \left[ V \right] – \left( {\frac{{\left[ r \right]^{4} }}{64} – \frac{{\left[ r \right]^{2} }}{4}\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right)\left\{ {\left[ C \right] + \left[ D \right]\left[ {B^{*} } \right]^{ – 1} \left[ {A^{*} } \right]} \right\}^{ – 1} \left[ {U_{0} } \right] \\ \end{aligned}$$

(47)

Having established \(\left[ \varphi \right]\), \(\left[ \mu \right]\) and \(\left[ Z \right]\), it is possible to formulate quintic spline equations \(S_{i} \left( x \right)\) for \(i = 1, 2, \cdots , N – 1\). \(S_{i} \left( x \right)\) 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}\), \(\eta_{i}\) for \(i = 1, 2, \cdots , N – 1\) must first be determined from (21), and based on the assumed boundary condition \(\eta_{N} = 0\). Having all \(\eta_{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} \left( x \right)\) for \(i = 1, 2, \cdots , N – 1\). The only difference is that, after determining \(\eta_{i}\) for \(i = 1, 2, \cdots , N – 1\), the value of \(\eta_{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,32,33}.