Calculation of the Absorbed Electron Energy 3D Distribution by the Monte Carlo Method, Presentation of the Proximity Function by Three Parameters α, β, η and Comparison with the Experiment

The paper presents a program for simulating electron scattering in layered materials ProxyFn. Calculations show that the absorbed energy density is three-dimensional, while the contribution of the forward-scattered electrons is better described by a power function rather than the commonly used Gaussian. It is shown that for the practical correction of the proximity effect, it is possible, nevertheless, to use the classical two-dimensional proximity function containing three parameters: α, β, η. A method for determining the parameters α, β, η from three-dimensional calculations based on MC simulation and development consideration is proposed. A good agreement of the obtained parameters and experimental data for various substrates and electron energies is shown. Thus, a method for calculating the parameters of the classical proximity function for arbitrary layered substrates based on the Monte Carlo simulation has been developed.


Introduction
One of the most common methods for creating micro-and nanostructures in microelectronics is electron-beam lithography (EBL). Although EBL is less productive than photolithography, it turns out to be very convenient to create small structures consisting of elements of very different sizes. This makes it in demand when creating micro-and nano-objects for scientific research. Such structures can be used for studies of superconductivity [1,2], X-ray radiation [3][4][5], electrophysical properties of various materials (for example, graphene [6,7]) and in many other areas of physics.
One of the features of EBL is the effect of electrons back-scattered in the substrate on the electronic resist. In this case, they carry out additional exposure of the electron resist on an area usually much larger than the size of the primary electron beam. This effect is commonly called the "proximity effect" [8]. Taking into account the influence and correction of the "proximity effect" on the dose absorbed by the electronic resist allows increasing the accuracy of electron lithography, shortens the time to fabricate structures and increases the yield of a suitable product, since it reduces the sensitivity of lithography to random errors. To correct the influence of the "proximity effect" when calculating the exposure dose, the proximity function (PF) is used, i.e., the distribution of the absorbed energy in the electron resist during electron-beam scattering.
The classical PF I (x,y) consists of two Gaussians [8], does not change over the depth of the resist film z, is dimensionless and is normalized to unity [9,10]. It is written in the following form: where r 2 = x 2 + y 2 , η is the ratio of the total energy left by the reflected electrons to the energy of the forward-scattered electrons. The first Gaussian in (1) describes the distribution of the energy left in the resist by the forward-scattered electrons (characterized by parameter α), and the second one describes the distribution of the energy left in the resist by the back-scattered electrons (characterized by parameter β) [11,12]. From (1), it follows that for the practical use of PF it is necessary to know the values of the parameters α, β, η. The experience of practical correction [10,13,14] and extensive simulations [15][16][17] show that the three parameters found in the experiment are in most cases quite enough to obtain the required lithography accuracy.
There is a large number of works devoted to calculations [12,[18][19][20] and to experimental measurements of PF [21][22][23]. The MC method [24] has long been used for determining the parameters of different functions describing the proximity effect [12,[18][19][20] by fitting the functions to the MC calculated spatial (3D) distribution of absorbed energy. Figure 1 shows the distribution of the absorbed energy density G (r, z) for a PMMA film of thickness H 0 = 1um on a silicon substrate, calculated by the Monte Carlo method. It can be seen that G (r, z) varies in the thickness of the resist and is really a 3D function. For fitting the 2D proximity function to 3D data, usually a distribution of absorbed energy at the boundary resist-substrate (i.e., G (r, z = H 0 )) was used [12,[18][19][20]. However, our experience has shown that consideration of the distribution only at the resist-substrate interface or averaging over the resist thickness does not allow one to obtain PF parameters that are in good agreement with the experimental values [21]. A particularly large discrepancy arises in the determination of the η parameter. Obviously, the cause is ignoring the process of development. Experimental methods for determining the PF parameters inevitably include development processes in the measurement procedure. When considering experimental methods, we prefer methods using resist development time ("vertical" methods) [21,22], over methods using measurements of the transverse dimensions of test features such as line widths or ring widths (e.g., [23]) ("horizontal" methods). Experimental determination of PF parameters is a laborious process and can take [21] several days or even weeks. The Monte Carlo calculation takes a few minutes.
Therefore, the purpose of this work will be to develop an exclusively computational method for obtaining the parameters α, β, η of the 2D classical proximity function from the 3D absorbed energy density calculated by the Monte Carlo method with careful consideration of development. For comparison, experimental data α e , β e , η e will be taken from [21].  Therefore, the purpose of this work will be to develop an exclusively computational method for obtaining the parameters α, β, η of the 2D classical proximity function from the 3D absorbed energy density calculated by the Monte Carlo method with careful consideration of development. For comparison, experimental data αe, βe, ηe will be taken from [21].

Calculation of Electron Scattering in Layered Materials by the Monte Carlo Method
We have developed an algorithm and implemented it in the ProxyFn program for the fast simulation of electron scattering in layered materials [25,26]. The structure of the algorithm is close to the one described in the work [24], but all expressions for the calculation are taken from the book by L. Reimer [27].
In the Monte Carlo simulation, it is assumed that electrons "move" in a straight line, continuously losing energy, until elastic scattering. A new direction of electron "moving" is played out using a screened Rutherford cross section [27]. The length of the straight segments is chosen randomly based on the total cross section of elastic scattering along the trajectory. Before scattering, in correspondence to the continuously slowing down approximation, the electron energy decreases, taking into account the length of the segment and the current value of stopping power (for details see [12][13][14][15][16][17]24]).
In calculations, the sample is a layered structure consisting of an arbitrary number of layers and arbitrary materials. For convenience, the sample is automatically divided into cells by a grid. In the cells of the partition grid (ri, zj), the energy of electrons E (ri, zj) left by them when passing through these cells and the number of stopped electrons N (ri, zj) are remembered. Additionally, the coefficients of reflection and transmission of electrons

Calculation of Electron Scattering in Layered Materials by the Monte Carlo Method
We have developed an algorithm and implemented it in the ProxyFn program for the fast simulation of electron scattering in layered materials [25,26]. The structure of the algorithm is close to the one described in the work [24], but all expressions for the calculation are taken from the book by L. Reimer [27].
In the Monte Carlo simulation, it is assumed that electrons "move" in a straight line, continuously losing energy, until elastic scattering. A new direction of electron "moving" is played out using a screened Rutherford cross section [27]. The length of the straight segments is chosen randomly based on the total cross section of elastic scattering along the trajectory. Before scattering, in correspondence to the continuously slowing down approximation, the electron energy decreases, taking into account the length of the segment and the current value of stopping power (for details see [12][13][14][15][16][17]24]).
In calculations, the sample is a layered structure consisting of an arbitrary number of layers and arbitrary materials. For convenience, the sample is automatically divided into cells by a grid. In the cells of the partition grid (r i , z j ), the energy of electrons E (r i , z j ) left by them when passing through these cells and the number of stopped electrons N (r i , z j ) are remembered. Additionally, the coefficients of reflection and transmission of electrons from the entire sample are determined, as well as the coefficients of the absorption of electrons and energy in all layers. To start the calculation, it is necessary to know only the starting energy of electrons, film thickness, chemical formulas of materials and their density. For many elements and materials, the chemical composition and density can be selected from the built-in database.
To speed up the calculation, a cylindrically symmetric and nonuniform grid r i , z j with center r = 0 on the beam axis is used. The z axis is perpendicular to the layers and directed from the outer boundary of the resist to the sample. z = 0, r = 0 is the point of entry of electrons into matter. The partition grid is set automatically.
The calculation does not consider the generation, scattering and absorption of secondary electrons, to which the fast electron gives up energy during deceleration. Although it is the secondary electrons with energies up to 50 eV that do all the "work" of breaking chemical bonds or forming bonds in resist molecules, their track length in materials does not exceed 10 nm [28] and can be taken into account by the convolution of the calculated absorbed energy density with the corresponding Gaussian. The Monte Carlo calculations also do not take into account the charging of dielectric layers, which can significantly change the trajectory of an electron. We believe that this problem can be effectively dealt with in electron lithography (for example, by applying a conductive film to the resist and grounding it, or by using a conductive resist). The PF calculation time for ten thousand trajectories takes several minutes.

Integral Proximity Function: Fitting of Absorbed Energy Distribution by Elementary Functions
It is not very difficult to implement the Monte Carlo algorithm for simulating electron trajectories. Difficulties arise when analyzing the calculation results and fitting the simulation results. After enumerating a large number of options, the following function was chosen to interpolate the density of the distribution of the electron absorbed energy G(x,y,z) in the sections z = const: where r is the distance to the beam axis. The first (delta-shaped) element in expression (2) describes the electrons of the primary beam, which have not experienced one scattering on atomic nuclei. The second and third elements in (2) describe (on a qualitative level) singly and multiply scattered electrons, respectively. The coefficient C δ (z) in the first approximation decreases exponentially with the penetration depth z, as where L f is the free length. It is inversely proportional to the total cross section of electron scattering in the resist and is equal to several tens of nanometers (about 80 nm for 25 keV electrons in PMMA). In the experiment, it is not easy to separate the first and second elements in expression (2) due to the fact that the initial electron beam is not delta-shaped. Note that the fitting parameters C δ (z), C a (z), C b (z), α(z) and β(z) depend essentially on the depth z and, in this case, the relation α(z) << β(z) is fulfilled.
To search for the fitting parameters C δ (z), C a (z), C b (z), α(z) and β(z), it is convenient to use not the distribution density G(x,y,z) itself, but the integral density of the absorbed energy G r (r,z) obtained from the expression: Note that the integral density G r (r,z) can be given a physical meaning. Consider a special structure in the form of an infinite plane with a cut out circle of radius r. It turns out that exposure of such a structure with a single dose leads to the absorbed energy density dE/dz at the center of the circle just equal to G r (r,z).
Using (2)-(3), we obtain for the three-dimensional proximity function: here On the other hand, the integrated absorbed energy density I r (x,y) in the case of the classical proximity function (1) will be equal to: Figures 1 and 2 show an example of calculating the integrated absorbed energy density G r (r,z) and its approximation for a 1 µm thick PMMA e-beam resist film on a 300 µm thick silicon substrate at an initial electron energy of 25 keV. The number of trajectories considered is 200,000. The classical PF consisting of two Gaussians does not approximate G r (r,z) very well, and function (4)

Fitting with Three Parameters: Analogue of Experiment
An experimental method for determining the classical PF I(x,y,α,β,η) (1) was proposed in [21]. The idea of the method is to search for such parameters (б, в, з) so that the test structure, consisting of elements of different sizes, after correcting the proximity effect (calculating the corrected dose based on the classical PF I(x,y,α,β,η)) and exposure the positive resist, will be revealed exactly to the substrate in the center of each element. This experimental method and the measured parameters б, в, з from [21] have been used for more than 20 years to correct the proximity effect in the NanoMaker software and hardware complex for electron-beam lithography (www.nanomaker.com, accessed on 26 May 2022) with consistently good results.
A similar method to search for the parameters of the classical PF I(x,y,α,β,η) is used in this work. As in the experimental method, the calculation of the exposure dose T(x,y) (electron density per unit area) is based on the classical PF I(x,y,α,β,η) from expression (1), but the simulation is performed instead of actual exposure and development. The ab- Thus, the depth-dependent three-dimensional PF G(r,z) can, in principle, be used to correct the proximity effect in e-beam lithography. However, a three-dimensional PF has significantly more fitting parameters than a classical PF with only three fitting parameters (б, в, з). This leads to a complication of calculations. On the other hand, as mentioned, the practical use (to correct the "proximity effect") of the classical PF leads to good results. Therefore, our next purpose is to determine the effective parameters of a two-dimensional PF from the three-dimensional Monte Carlo simulation results.

Fitting with Three Parameters: Analogue of Experiment
An experimental method for determining the classical PF I(x,y,α,β,η) (1) was proposed in [21]. The idea of the method is to search for such parameters (б, в, з) so that the test structure, consisting of elements of different sizes, after correcting the proximity effect (calculating the corrected dose based on the classical PF I(x, y, α, β, η)) and exposure the positive resist, will be revealed exactly to the substrate in the center of each element. This experimental method and the measured parameters б, в, з from [21] have been used for more than 20 years to correct the proximity effect in the NanoMaker software and hardware complex for electron-beam lithography (www.nanomaker.com, accessed on 26 May 2022) with consistently good results.
A similar method to search for the parameters of the classical PF I(x,y,α,β,η) is used in this work. As in the experimental method, the calculation of the exposure dose T(x,y) (electron density per unit area) is based on the classical PF I(x,y,α,β,η) from expression (1), but the simulation is performed instead of actual exposure and development. The absorbed dose (density of absorbed energy per unit volume) D(x,y,z) in the simulation is calculated based on the three-dimensional PF G(x,y,z) obtained by Monte Carlo from expression (2). For the dimensionless classical PF I(x,y,α,β,η), the distribution of the absorbed dose D(x,y) is as follows: The development of a positive e-beam resist is simulated in the approximation of isotropic, local etching [29,30]. Then, the development rate V can be written as follows: where γ is the contrast of the resist; and D 0 , V 0 are the technological constants. For a positive e-beam resist, the sensitivity T 0 is defined as the exposure dose at which an element with dimensions much larger than в is revealed in the center exactly to the substrate. A brief description of the approach presented in Appendix A is as follows. The proposed method consists in considering a number of circular elements of different sizes R. Exposing a circle with a uniform exposure dose results in an absorbed dose distribution with a maximum exactly at the center of any circular element at all resist depths z. From isotropic local etching theory [26,30], it follows that the development front reaches the substrate for the first time namely at the center of the round. Development times T R and T i R calculated for two different exposure models (for classical PF (5) and for three-dimensional PF (4)) are dependent on element radius R. Due to the simplicity of the elements, these times can be easily calculated by formulas. Further, such parameters of the classical PF б, в, з are searched with a special procedure that minimizes the objective function (9) using the ratio of the times T R /T i R . To search for б, в, з, a set of 10 round elements of radius R n (n = 1, . . . ,10) was used. The R n value varied from the minimum value α(z) to the maximum value β(z), 0 ≤ z ≤ H 0 , where α(z), β(z) are the interpolation parameters of the integrated density of the absorbed energy G r (r,z).
The Appendix (A4) shows that the exposure dose ratio does not contain technological parameters D 0 , V 0, t 0 and can be calculated relatively easily. The exposure dose ratios (8) were calculated (for the given parameters б, в, з) for the entire set of circles R n , and the objective function was composed from them.
Further, the values of the parameters α s , β s , η s were determined by minimizing S(α, β, η). This method allows us to calculate quickly the parameters of a classical PF, depending on the thickness of the resist, the type of substrate (including the possible layered structure of the substrate), electron energy, etc.

Results and Discussion
A comparison of the experimental parameters α e , β e , η e and those obtained from the simulation based on Monte Carlo calculations α s , β s , η s бs is necessary to verify the correctness of our chosen physical and mathematical models of electron scattering and resist development, as well as to verify the accuracy of calculations. It seems to us that a quantitative assessment of the accuracy of the Monte Carlo calculation (at least for thick resists) has not been performed before.
The experimental data were taken from [21]. In the experimental method used in this work, each of the three parameters α e , β e , η e was measured in a separate test. All tests used a PMMA electronic positive resist (chemical formula CH 2 C(CH 3 )(COOCH 3 ), density 1190 kg/m 3 ). In calculations based on the Monte Carlo method, all three parameters α s , β s , η s were searched simultaneously.
First, the data on the beam size б will be compared, and then the results of calculating в and з. The value of б depends on the energy, thickness and material of the resist and does not depend at all on the type of a substrate. In addition, the experimental value of б e is influenced by the initial beam size б 0 , which is determined by the electron microscope setting (focusing, astigmatism) and beam jitter. In fact, the value of б s obtained by the Monte Carlo method should be compared with parameter α 2 e − α 2 0 . Table 1 shows a comparison of the size of the forward-scattered electron beam б, obtained from the б e experiment (except the initial beam size б 0 ) and calculated on the basis of the Monte Carlo method б s for three energies E = 15, 25 and 35 keV and a set of resist thicknesses H 0 = 100, 200, 500, 1000 and 1500 nm on a silicon substrate. The experimental data б e were interpolated by the formula α 2 e = A E H 3 0 /E 2 + α 2 0 , where the constants A E and the initial beam dimensions б 0 for different energies E were obtained from the experiment [21]. In the Monte Carlo calculation, the beam was assumed to be absolutely thin. Table 1 shows that the experimental data for б e and the results of calculating б s are in good agreement. Table 1. Comparison of the parameters of the proximity function αs, (obtained by the method described above based on the Monte Carlo calculation) and αe, (calculated from the results of interpolation of experimental data [21]), for three values of the electron energy E and different thicknesses H 0 of the PMMA resist. The substrate is Si. A comparison of the calculated (β s , η s ) and experimental (β e , η e ) for various substrates and accelerating voltages can be carried out using the data in Tables 2-4. For all cases, a positive PMMA resist 500 nm thick with a contrast γ = 3 was used. The values of α s , β s , η s turned out to be stable with respect to the change in contrast, and hardly changed for г = 2.5, 3 or 4. For comparison, the experimental data for Si, GaAs, Al 2 O 3 and mica from [21] were used; the data for Ge and C (diamond) substrates were specially measured for this work by the method [21]. The calculated values of β s with an accuracy ofŷ10% coincided with the experimental values of β e , for з the accuracy wasŷ25%.   Note that the parameters в and з did not depend on focusing (if α s << β s ) and were determined only by the properties of the resist, substrate and the initial energy of electrons; therefore, they have a fundamental value.

Conclusions
The algorithm for the fast calculation of the absorbed energy density of electrons G(r) by the Monte Carlo method for layered materials was described.
To interpolate the calculated absorbed energy density of electrons G(r), a fitting function (2) was proposed, which described well the distribution of the absorbed energy of electrons in layered materials depending on the distance to the center of the beam r and on the depth z. The power member describing the scattering of primary electrons seemed to be nontrivial and was considered for the first time.
A numerical procedure was proposed that takes into account the development of the resist and makes it possible to replace the complex 3D distribution of the absorbed energy with a classical (two-dimensional) proximity function with three parameters б, в, з.
The examples of calculating the parameters α s , β s , η s of the proximity function were shown for different energies of electrons and substrates and their comparison with the experimental data α e , β e , η e . Calculations of в with an accuracy ofŷ10% coincided with the experiment; for з, the accuracy wasŷ25%.
Thus, it can be argued that the experimental confirmation of the accuracy of calculating PF by the Monte Carlo method and the procedure of interpolating PF with three parameters (б, в, з) has been obtained.  Data Availability Statement: Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The ratio of exposure doses is found for classical and three-dimensional proximity functions. The density of absorbed energy (dose) in the center of a circle with the radius R, exposed with a constant dose T (density of electrons per unit area), is expressed through the integral proximity function: D(R, z) = (G r (0, z) − G r (R, z)) T (A1) For the classical dimensionless PF I(x,y,б,в,з), the absorbed energy distribution is written as follows: = (I r (0, α, β, η) − I r (R, α, β, η)) T T r = (1 − I r (R, α, β, η)) T T 0 (A2) In order for the resist to reveal exactly in the center of the circle to the substrate, the following condition should be held: D/D 0 = 1.
Then, it follows from (A2) that the exposure dose T R should be equal to: The "ideal" exposure dose T i R is calculated using the three-dimensional PF. In the center of each element, the absorbed dose has a maximum on the XY plane in each section z = const; therefore, a positive electron resist is revealed vertically along the z axis in the center of the circle [30]. The absorbed dose D R (z) in the center of a circle with the radius R n exposed with a certain dose T is obtained from expression (A1). The time t for the development of a positive resist with a thickness H 0 to the substrate in the center of the circle is t = H 0 0 dz V(z) . Then for a circle with the radius R: The "ideal" exposure dose T i R can be obtained from expression (A3), at which a circle with the radius R in the center is revealed to the bottom for a given development time t 0 : The sensitivity of the positive resist T 0 can also be calculated using this expression at R = ∞: As a result, the ratio of exposure doses T R /T i R is obtained, which does not contain unknown technological parameters D 0 , V 0 and t 0 :