## 1. Introduction

Increasing the electronic circuit density is a major trend in microelectronics. Lithography is the key technology for the fabrication of very large integrated circuits with smaller device sizes. Today conventional optical lithography is approaching its fundamental physical limits. Electron and ion lithographies are among various candidates capable of sub-150 nm resolution for the new generation of lithographic techniques. Nanometer scale device fabrication rules require tight control of the developed polymer resist profile. Process simulation is a key tool for optimization of the obtained lithography results.

The goal of the computer simulation of the processes at electron and ion beam lithography (EBL and IBL) is the resist profile prediction of developed patterns after exposure of samples, covered by a sensitive polymer resist layer, which is sensitive to irradiation by accelerated particles. The accuracy of the simulated resist profiles strongly depend on the physical knowledge of the processes as well as on the accuracy of the process parameters.

The main step of such modeling:

## 2. Deposited energy in the case of electron beam lithography simulation

The first main step of a complete mathematical model for electron beam lithography (EBL) simulation is to simulate the exposure of polymer resist films. During the exposure process the resist material modifies the local solubility rate. The main goal during the exposure process modeling is the calculation of the absorbed energy space distribution. There exist two types of method to calculate the deposited energy in the resist (latent image): analytical methods (Hatzakis et al., 1974; Hawryluk et al., 1974; Raptis et al., 1998) and numerical methods (Kyser & Murata, 1974; Adesida et al., 1979; Vutova & Mladenov, 1994). The analytical method is based on some particular approximations that simplify the nature of the real process (small-angle and diffusion particle scattering, single-component targets, point source or source of homogeneous cross section, etc). These assumptions do not correspond to the real process of the beam scattering within the target. Nevertheless, the simulation of energy deposition in the resist film coated on patterned substrates using analytical methods is very difficult if not impossible.

The numerical methods, based on the Monte Carlo (MC) technique for the statistical electron trajectory modeling and energy loss calculation, have been extensively developed and have become the most accepted methods in this field. The MC method mirrors the real process and in the case of large-number trajectory modeling, assures high statistical accuracy and satisfactory consistency. On the other hand, the MC method is ideal for parallel processing computers.

In the proposed simulation model the absorbed energy space distributions are calculated using MC algorithm for electron penetration and energy-loss calculation, which has four sub-steps:

forming an electron scattering model and calculating the discrete absorbed energy distribution in the resist film due to a point beam,

approximating the absorbed energy using analytical expressions (energy deposition function EDF(r,z)),

convoluting the function EDF(r,z) with the actual current distribution in the electron beam used for exposure,

calculating the spatial distribution of the absorbed energy density in the resist which determines the obtained latent image during the electron beam exposure process of the desired layout.

The process resolution is limited by the phenomena of forward deflection and backscattering of the electrons during their passage through the resist layer and the substrate. When a real micro-image is exposed, the absorbed energy in every resist volume point can be calculated by summing up the energy losses, obtained in regions far from the beam incident point (several μm). In this way the scattering of the beam electrons limits the lateral resolution of the exposed lithographic patterns and patterning of dense high-resolution layouts. This phenomenon is known as a proximity effect and the function EDF(r,z) characterizes the so-called proximity effect (undesired exposure dose due to backscattered electrons) (Chang, 1975).

Calculation of energy deposition in the resist film due to a point beam (with a zero-width beam diameter δ-function) requires to investigate thousands electron trajectories, i.e. millions collisions between accelerated electrons and scattering target atoms (an elementary collision sequence with target atoms). We calculate the characteristic changes in the particle motion for each collision, assuming a straight-line trajectory between two collisions (Fig.1).

The scattering atom is presented by a shielded Coulomb potential. The scattering angle θ of the penetrating electron is calculated using the differential scattering cross-section for the penetrating electron (for elastic collisions with target atoms) and assuming a Rutherford shielding potential presenting the scattering atom:

where σ is the differential scattering cross section [mm^{2}], Ω is the solid angle, V is the velocity and E [eV] is the energy of the penetrating electron, Z is the atomic number of the scattering atom. The shielding parameter β characterizes the minimal scattering angle at which the value of dσ/dΩ does not increase more when the value of θ decreases:

and then the total cross section of the elastic scattering is:

where σ [m^{2}] if E [eV].

The non-elastic scattering weakly changes the electron trajectory due to its domination only at very small values of the scattering angles. It is taken into account as Z^{2} is replaced by Z(Z+1) in Eq.(1) and Eq.(3). The mean energy losses of the penetrating electrons, evaluated per one unit of the trajectory length are:

where dE/dξ [eV/Å], the target density ρ [g/cm^{3}], E [eV], N_{A} is the Avogadro’s number, c_{i} is the weight portion of the i-kind of the target atoms, M_{i} is their atomic weight, I_{m} is the average ionization potential.

The program involves a MC technique to calculate the kind of the atom taking part in the collision, the azimuthal angle value, ect. The probability for a collision with the target atom of k-kind is:

where m is the number of the different kinds of the target atoms. The scattering cross-section σ_{i} is calculated using Eq.(3) for the electron energy E before the collision. The concentration n_{k} of this kind of atoms is calculated using:

where

where

where R_{1} is a random number evenly distributed in [0,1]. The energy losses of the penetrating electron at the interaction with the target electrons along this path are calculated using the Bethe energy-loss Eq.(4). The scattering angle of the collision is calculated:

The azimuthal angle φ is given by:

where R_{2} and R_{3} are also random numbers. Due to the big difference concerning the weights of the particles, the energy losses at the collision between the penetrating electron and scattering atom are neglected. The values of *Δξ, θ* and *φ* for each collision are calculated. The calculation is repeated to yield new position in the target for a new set of random numbers until the energy of the electron falls below a predefined value (500 eV) or the electron leaves the target. Then the electron trajectory is calculated (Fig.2).

The trajectory can be presented in the coordinate system Oxyt_{z} (easily connected with a cylindrical system, and *r=(x*
^{
2
}
*+y*
^{
2
}
*)*
^{
1/2
}), where the axis *t*
_{
z
} is parallel to the initial direction of the penetrating electron motion (usually this direction is perpendicular to the target surface) and using the relation:

The angles *ψ*
_{
i-1
}
*, ψ*
_{
i
} and *θ*
_{
i
} are shown in Fig.1. The angle *φ*
_{
i
} is the azimuthal angle of scattering for the i-th collision. The point depth t_{z} for the (i+1)-th collision is calculated using the formula:

The distance to the coordinate axis *0t*
_{
z
} is calculated by analogy assuming that the angle *χ*
_{
i
} determines the trajectory path projection on the axis *0r,* as well as that the angle φ_{
i
} corresponds to the azimuthal angle *φ*
_{
i
}. The following equations present the relations between these angles:

Then:

The energy losses for each segment (a straight line trajectory between two collisions) of the free path *Δξ*
_{
i
}, is represented as

and can be assumed as the energy losses in the point of the i-th collision for the corresponding elemental volume in the target. We check the energy of the electron and its position in the target. The calculation for this particle stops if the energy falls below the minimal predefined value (500 eV) or the electron leaves the sample. Then the calculation for new electron starts. After summing up the losses corresponding to the all *N* penetrating electrons, we obtain discrete data (a two-dimensional data array) for a radial energy deposition function (EDF) at various resist depths from the point of beam incidence. It is possible to use various cell dimensions at different radial distances: lower values near the point of beam incidence and higher values far away from this point. When using a Cartesian coordinate system instead of a cylindrical one it is possible to simulate beam incidence inclined to the resist surface (Gueorguiev, 1996). There are procedures for re-calculation of the free path and electron stopping power when a penetrating electron crosses the interfaces in multilayer structures (Gueorguiev et al., 1994). To achieve a satisfactory statistical accuracy, using MC calculations, a large number of particle trajectories are simulated (10-20 x 10^{3}). A detailed description and obtained results are presented in (Vutova & Mladenov, 1994; Mladenov & Vutova, 2002; Vutova, 2007).

### 2.1. Results and discussion

We have developed a program which realizes the MC algorithm to model the electron scattering in multilevel multicomponent amorphous targets, named TREM-MV. Using the simulation tool, calculations over a wide range of primary electron beam energies and resist thickness values are performed for different substrates (Vutova & Mladenov, 1994; Gueorguiev et al., 1995; Mladenov & Vutova, 2002; Vutova et al., 2007). In Fig.2 calculated trajectories in the case of a 1 μm thick poly-methyl methacrylate (PMMA) (a widely used polymer resist) on Si substrate are shown. The chemical composition of this resist is C_{5}H_{8}O_{2}, the efficient atomic number is *Z*=3.6, the atomic mass is *M* = 6.7 and the polymer density is ρ= 1.22 g/cm^{3}.

Fig.3 presents the radial distribution of the electron energy deposition in the resist film for 0.12 μm thick PMMA on Si due to a point beam by tracing 10 000 electron trajectories (for each simulation). Two characteristic regions can be seen: a narrow one, near the beam axis, with an abrupt drop, which represents the forward scattering electron contribution to the function EDF(r,z) and an wide part with a slightly varying drop, which corresponds to the backscattering electron contribution. As the initial electron energy increases, the area of the forward scattering electrons broadens insignificantly, while the region of the backward scattering electrons undergoes considerable broadening.

Another presentation of the calculation results is by means of the equi-energy density contours as shown in Fig.4. When increasing the polymer thickness (Fig. 4(a)-(c)) a deviation of the most distant equi-energy density contours from the beam axis is observed, which is due to broadening of the forward scattering electron area. Further increase of the resist thickness (Fig.4(d)) results in a deviation of the equi-energy density contours back to the beam axis which is due to the beam intensity decrease. The doses corresponding to the presented equi-energy density contours decrease starting from the beam axis as follows: 1; *1/2*; *1/4; 1/10; 1/20; 1/75.* These results can be used both for the latent image profile prediction under particular exposure conditions and for the exposure condition optimization, allowing the desired profile to be obtained. For example, one can see from Fig.4, that to obtain a latent image with vertical walls (primary e-beam energy, 20 keV) it is preferable to use thin polymer film whose thickness is less than 0.4 μm.

The high resolution of EBL however may be degraded by the lateral scattering of incident electrons which causes undesired exposure of unintended regions of the resist. This phenomenon is commonly referred to as proximity effect. In the case of patterning the thin films of the most widely used high temperature superconducting (HTS) material, namely YBa_{2}Cu_{3}O_{7-δ} (YBCO) deposited on SrTiO_{3} (STO), MgO, ZrO_{2}: Y_{2}O_{3}, LaAlO_{3}, NdGaO_{3}, etc. substrates an enhanced proximity effect has to be taken into account because of their

retively high effective atomic numbers (Gueorguiev, 1994). Using our program for MC simulation of the processes of penetration and scattering of accelerated electrons in solids, the radial distributions of the EDF in 125 nm PMMA resist layer coated on structures YBCO thin film/substrate are obtained for an e-beam in the form of a zero-widt δ function, 30 000 electron trajectories, and the following parameters:

the substrate material (STO and MgO),

the e-beam energy

*E*_{ 0 }(25, 50 and 75 keV), andthe YBCO film thickness

*d*(0, 100, and 300 nm).

In Fig.5, the calculated radial distributions of the absorbed electron energy density at *E*
_{
0
} = 75 keV are shown as an example. The results show that the EBL on the above mentioned targets is associated with an enhanced proximity effect in comparison with that on the conventional in microelectronics targets PMMA/Si substrate or PMMA/SiO/Si substrate. Moreover, the HTS thin film causes an additional backscattering of penetrating electrons and, hence, an additional proximity effect (in comparison with the targets PMMA/STO and PMMA/MgO) in the regions close to the incident point of the electron beam. This effect is as greater, as thicker is the film, as lighter is the substrate, as lower is the beam energy and is not completely eliminated even at energies as high as 75 keV, especially for the film thickness 300 nm, as well as for the lighter substrate (MgO).

## 3. Approximation of the discrete data for the deposited energy in the resist using analytical functions

Due to large lateral scattering of the penetrating electrons, the exposure of many image segments effects the total deposited energy in a specific resist point. This effect, known as a proximity effect, requires high accuracy evaluation of the EDF at large radial distances i.e. far away from the point of incidence. It should be pointed out that the increase of the number of trajectories being modeled, with the purpose of achieving statistical consistency for large lateral distances (characterizing the backward scattered electrons and giving the greatest contribution to the proximity effect) is not quite effective, as only a few trajectories travel through these regions. A MC methodology and a corresponding computer program BET-MK are developed for transformation of the numerical data array, representing the absorbed energy space distribution when exposing one point from the resist surface, into the form of analytical functions (Vutova & Mladenov, 1994). In this way the problem concerning the insufficient statistics of the discrete data for the absorbed energy in the case of large lateral distances is overcome. The main features of this methodology are as follows:

Transformation of numerical data for the absorbed energy discrete space distribution in the form of analytical functions. The data arrays containing the energy distribution when exposing one point from the resist surface are obtained using the TREM-MV program.

The absorbed energy at some resist depth in the case of Si substrate is approximated as a sum of two Gaussians:

where

*k*is a normalization constant,*δ*_{ f }and*δ*_{ b }are the characteristic widths of the forward and the backward scattering particles, and*η*_{ E }is the ratio of the energy depth dissipation of the backward scattering particles to that of the forward scattering particles. The input data for the program BET-MK is the 2-D arrays EDF(r,z) containing the absorbed energy distribution values obtained as a result of the trajectory modeling. The first Gaussian (with standard deviation*δ*_{ f }) dominates for the short lateral distances and describes the energy deposition from the forward scattered electrons. The second Gaussian (with standard deviation*δ*_{ b }) dominates for the long lateral distances and describes the contribution from the backscattered electrons. The parameters*δ*_{ f }*, δ*_{ b }and*η*_{ E }are called proximity effect parameters.The parameter values

*(δ*_{ f }*, δ*_{ b }*, η*_{ E }*)*are calculated using an original MC technique, instead of the non-linear least-square method and an arbitrary kind of distribution. The technique comprises the mean square deviation minimization by the interval length decrement for each of the parameters chosen. The minimization is made in an iteration loop.

The main advantages of the MC technique described above are:

in contrast to some of the least-squares methods, it does not allow the possibility of an infinite loop in the case of a local minimum;

it enables to approximate an arbitrary kind of distribution of numerical data with a corresponding analytical function.

Fig.6 presents the EDF in the resist film for 0.4 μm thick PMMA on Si for 50 keV beam energy. The standard deviation for forward *(δ*
_{
f
}
*)* and backscattering *(δ*
_{
b
}
*)* contributions (Eq.19) as calculated from the simulation data are: *δ*
_{
f
}
*=* 0.165 μm*, δ*
_{
b
}
*=* 7.769 μm*, η*
_{
E
} = 1.2_{.}

The values of the parameters *δ*
_{
f
} and *δ*
_{
b
} for various electron beam energies and different resist (PMMA) thicknesses are shown in Figs.7-8. In the triangular diagram the points corresponding to the values of *δ*
_{
f
} or *δ*
_{
b
}, the resist depth, and the beam energy represent vertices of inscribed triangles. Using these results one can approxima determine the *δ*
_{
f
} and *δ*
_{
b
} values for different beam energies and resist thicknesses. For instance, if the beam energy is 15 keV and the PMMA thickness is 0.8 μm*,* the *δ*
_{
f
} value at the depth of 0.1 μm is in the range of [0.21, 0.29] μm*,* and that at the resist-substrate interface is between 0.24 and 0.33 μm*.* If the beam energy is 20 keV, the PMMA thickness is 0.7 μm*,* the *δ*
_{
b
} value at the depth of 0.1 μm is in the range of [0.13, 0.21] μm*,* and at the resist-substrate interface its value is between 0.14 and 0.23 μm*.* Similarly, if *δ*
_{
f
} (or *δ*
_{
b
}) values and the beam energy are known, one can determine the resist thickness and vice versa from known *δ*
_{
f
} (or *δ*
_{
b
}) and the resist thickness, the beam energy can be determined.

According to the simulation strategy, the next step in the modeling is the convolution of EDF(r,z) with the actual beam current distribution j(r,z), giving the single-energy deposition function (SEDF(r,z)). Let the electron beam is Gaussian distributed with a characteristic width δ^{*}
*,* i.e.

where *I*
_{
b
} is the incident beam current, *j* and *j*
_{
o
} are the current densities at a distance *r* and at the beam centre, respectively. Then the function for a real e-beam can be calculated by the convolution of (19) and (20). The result from this convolution is the function f^{*}
*(r),* that has an analytical representation in the form of (19), but its characteristic widths are now modified:

## 4. A three-dimensional model for absorbed energy calculation

To obtain the absorbed energy space distribution when exposing an arbitrary pattern, using an arbitrary exposure dose distribution, one must take into account the influence of a large number of exposed points. This is due to the fact that the lithography micro-patterning includes many irradiated points. Integrating the data obtained by the computer simulation of the exposure of each point from the resist surface, the absorbed energy space distribution in the case of an arbitrary pattern can be calculated. Due to the large number of calculations, a simplified procedure should be used to calculate the integral space distribution of the absorbed energy. The main features of the procedure proposed in (Vutova & Mladenov, 1991) are as follows.

The two-dimensional data array containing the absorbed energy values at some resist depth is presented as (19). If the electron (or ion) exposure is uniformly distributed over an area

*A,*then the energy density can be expressed as:If the area

*A*is a simple pattern (i.e. a line or a rectangle) the integral (23) can be calculated using the tabulated error function:In the case of a more complex pattern, it should be divided into simple parts and then the corresponding values of the absorbed energy should be subsequently summed up. The formulas obtained for the absorbed energy density, when exposing either a line, a line segment, or a rectangle (Vutova & Mladenov, 1991), are given below. These simple patterns are sufficient to compose an arbitrary figure.

The procedure takes into account the radial variation of the absorbed energy as well as its modification versus the depth of the resist. To calculate the

*δ*_{ f }*, δ*_{ b }*,*an*d η*_{ E }parameter values, the linear approximation along the resist depth is used. In the electron exposure case, two resist depths are used, namely the resist surface and the resist-substrate interface (Table 1). The proposed methodology is realized in a computer program.

Let we have Cartesian coordinate system *Oxyz,* such that the axis *Ox* and *Oy* are in the plane of the resist surface, the axis Oz is in its depth, and *x*
^{
2
}
*+ y*
^{
2
}
*= r*
^{
2
}. Let we expose a line segment and the axis *Oy* coincides with the line segment, whose end points are labeled *a* and *b* (Fig.9).

We are interested in the cross section that is perpendicular to the segment line and whose analytical expression is *y* = *c.* From (19) and (23) we obtain:

Similarly, we can obtain the expressions for the density of the absorbed energy in a point lying at the *y* = c line, when exposing a rectangle ABCD (Fig.10):

(30) |

where *erf(t,σ)* is calculated using (24).

The absorbed energy density distribution in the case of a more complex topological structure, where the proximity effect cannot be ignored is shown in Fig.11. The results are presented for three inter-elemental spacings of d=0.5, 1.0, and 1.5 μm. This simulation result clarifies the effects of an adjacent element and its configuration on the energy distribution.

## 5. Absorbed energy approximation in the case of multilayers samples and heavy substrates

The proximity effect in the case of patterning the thin films of the most widely used HTS material, namely YBCO, deposited on two typical substrates (STO and MgO) is investigated (Gueorguiev et al., 1995; Gueorguiev et al., 1996; Gueorguiev et al., 1998; Olziersky et al., 2004; Vutova, 2007). HTS samples represent a more difficult case study since the substrate consists of bulk substrate STO or MgO and a very thin YBCO layer on top (multilayer substrate). The existence of the thin YBCO film between the bulk substrate and the resist changes the scattering phenomena and has to be carefully taken into account. This effect becomes more important as film thickness increases.

For substrates with larger mean atomic number and density the addition of one more function is necessary. In the case of YBCO films over STO or MgO substrates the addition of an exponential function in equation (19) is found to be adequate (Gueorguiev et al., 1996). Thus EDF(r) could be approximated as:

where the third term describes the energy deposition in the mid-lateral distances. The values of the parameters of this function (called proximity function) *α, β, γ, η, ν* and *k* are calculated using the MC technique (described above 3.).

### 5.1. Results and discussion

The variables studied in our work are the substrate material (STO and MO), the initial energy of accelerated electrons *E*
_{
0
} (25, 50 and 75 keV) and the HTS film thickness *d* (0, 100, 200, 300 and 1000 nm). The values of the proximity effect parameters are evaluated from the fitting of EDF(r) with a sum of suitable functions (Eq.29) and their dependence on all investigated variables are discussed (Gueorguiev et al., 1995; Gueorguiev et al., 1996; Gueorguiev et al., 1998; Olziersky et al., 2004). The absorbed energy distributions obtained and the calculated parameters of the proximity function can be used in a proper proximity effect correction algorithm as well as in a resist profile development model. The proximity effects are usually compensated for by applying proper correction algorithms which adjust the exposure dose and/or the shape and size of the exposed pattern. For the realization of such algorithms precise data are required about the spatial distribution of absorbed electron energy density in the resist. This distribution quantitatively describes the proximity effects.

In Fig.12 a comparison is made between the radial distributions of the absorbed energy density obtained by MC simulation for the structure 125 nm PMMA resist film /300 nm YBCO HTS film /MgO substrate at three beam energies -25, 50, and 75 keV, and the corresponding analytical fits. It is well seen that the combination of double Gaussian and exponential functions is a good approximation of these distributions. Although not shown here, the double Gaussian (Eq.19), as well as the triple Gaussian were also tested but they were found to be not adequate, especially in the intermediate regions (Fig.12).

Figure 13 shows the analytical fit to the radial distributions of absorbed energy density obtained by MC simulation for the structures 125 nm PMMA resist film /0, 100, or 300 nm YBCO HTS films / STO or MgO substrates at 25, 50, and 75 keV, respectively. Since the aim is to investigate the proximity effects caused by YBCO film as well as by the substrate (STO, MgO) the backscattered exposure is of primary importance. For this reason here, in contrast to the Fig.12, a linear scale for the *x* axis is applied which although it compresses data points close to the beam axis (associated with the forward scattered electrons), it ensures a better distinction between the distributions in their intermediate and distant regions (associated with backscattered electrons).

The peaks of the distributions of the absorbed energy density are commonly attributed to the forward scattering of electrons or, in other words, to the single scattering of primary electrons into small angles in the resist. This scattering depends on the beam energy as well as on the material and thickness of the resist. In Figs.12-13 it is seen that the maximum values as well as the widths of the peaks decrease with increasing beam energy. This can be explained by both the more efficient scattering of primary electrons and the higher energy loss in the resist at lower energies.

The calculated values of the parameters of the analytical function and their dependencies on the beam energy and on HTS film thickness are presented in the form of triangular diagrams as well as of 3D diagrams that can be used for easy approximate determination of the parameters at different beam energies and YBCO film thicknesses – Fig.14 and Fig.15.

The results show that the additional backscattering of primary electrons and, hence, the proximity effect, caused by the HTS film in the regions close to the incident point of the electron beam are not completely eliminated even at energies as high as 75 keV especially for the film thickness 300 nm as well as for the lighter substrate (MgO). The HTS thin film reduces the backscattering from the underlying substrate and this reduction is as greater as thicker is the film as well as lower is the beam energy.

In Fig.16, EDF(r) simulation results with the MC method in the case of YBCO/MgO substrates are presented. It is obvious that as d_{YBCO} increases, the overall energy deposition approaches the bulk YBCO case. When the YBCO layer is thin, due to the its relatively higher scattering parameters (density, mean atomic weight, mean atomic number) in comparison with those for the MgO substrate, the EDF(r) extends to regions far away from the point of incidence. On the other hand, when the YBCO layer is thick, the backscattered electrons come mainly from this layer rather than the MgO substrate. The external proximity effect obtained in regions far from the point of beam incidence (more than 4–5 μm) increases

when the d_{YBCO}< T_{crit}, while the proximity effect is lower in regions near to the point of beam incidence (T_{crit} is the YBCO thickness that causes energy deposition equal to the EDF(r) for bulk YBCO). The exponential function dominates in the intersection area between the two Gaussian functions which in the 50 keV case occurs in the ~ 40 nm to 2 μm range.

## 6. Comparison of simulated results with experimental data

The fundamental quantity in the simulation of EBL is the energy deposition due to point beam exposure (EDF(r)) or a single pixel exposure (SEDF(r)). When this quantity is in good agreement with the corresponding experimental data then it is possible to obtain good simulation results for the final resist profiles provided accurate dissolution rates are known.

Two different modules for the energy deposition calculation EDF(r,z) are applied in the case of Si and HTS films. The first module is based on the MC algorithm (Vutova & Mladenov, 1994; Gueorguiev et al., 1998) and the second based on an analytical solution of the Boltzmann transport equation (AS) (Raptis et al., 1998). Comparison between these approaches for various energies and substrates was performed in order to test the adequacy of calculated set of data. In addition, results from both modules were compared with experimental data for thin resist films on Si and HTS substrates (Gueorguiev et al., 1995; Olziersky et al., 2004; Vutova et al., 2007).

In Fig.17 the SEDF(r) in the case of 400nm thick epoxy resist on Si substrate is presented. The symbols represent the experimental data while the two lines represent the simulation results for MC and AS methods after convolution with the experimental beam diameter. Due to the 50 keV e-beam energy and the Si substrate used in this comparison, the backscattering contribution is spread out over a long distance from the point of incidence. It is obvious that both simulation methods provide results with good agreement with experimental data in the whole (forward scattered and backscattered) range.

In Table 2 the calculated β values in all cases (energy, substrate composition) is presented using both MC and AS modules, and the agreement between the calculated results is good.

In Fig.18 the simulation data (MC and AS) are presented along with the experimental SEDF(r) for STO and YBCO/STO substrates. The simulation results were prepared after convolution of the simulated EDF(r) (at the resist/substrate interface) from a point beam with a Gaussian beam diameter 30 nm. The SEDF(r) is obtained from exposure with 50 keV electrons in the case of a 0.4 μm thick resist film. In all cases, the agreement between experiment and simulation is satisfactory. Even in the second case, where the substrate is not homogeneous but consists of two layers (a thin YBCO layer over bulk STO) the discrepancy is limited.

## 7. Exposure process simulation at ion beam lithography (IBL)

A new look at previously obtained experimental results and simulation models applicable in IBL is a necessary task due to the progress in the development of IBL tools (Kaesmaier & Löschner, 2000) and their applications concerning exposure of one-component resist PMMA (Miller, 1989; Mladenov et al., 2001) as well as multi-component resists including chemically amplified resists (CARs) (Hirscher et al., 2001). The simulation tools for sub-150 nm ion beam lithography require more adequate models of the ion scattering, the energy losses and the development process. Accumulation and analysis of calculated results and physical parameters included in the models permit deeper understanding of the processes, comparison with experimental data and optimization of the simulation tools. Ion lithography of PMMA of thickness up to 400 nm is studied.

### 7.1. Model and computer program TRIM-MV for simulation of ion scattering and exposure intensity calculation in the case of ion bombardment

Stopping powers and trajectories of the penetrating atom particles in the materials under ion bombardment are quite well understood, due to the progress of ion implantation, ion sputtering and secondary ion mass spectrometry. Computer programs such as MARLOW, TRIM, PIBER, SASAMAL (Robinson & Torrens, 1974; Biersack & Haggmark, 1980; Adesida et al., 1982) were developed between 1970 and 1980. The program structure is the same as for electrons. The major difference is that the energy transfer is executed by two types of losses: nuclear stopping power, due to penetrating particle collisions with the sample atoms and electronic stopping power, due to penetrating projectile collisions with the electrons in the bombarded material. The first applications of these programs have been for PMMA (Adesida et al., 1982; Mladenov, 1985; Mladenov et al., 1985) to evaluate the ranges and energy losses.

The main computer program for ion scattering in amorphous materials using the Monte Carlo calculations is TRIM, first published by Biersack and Haggmark (Biersack & Haggmark, 1980), but later extended fully to cover collision cascades including the motion of recoils (Biersack & Eckstein, 1984). For the particular energy range in ion lithography the calculations are further speeded up by introducing an approximate but rather accurate treatment of the individual collision process. The program has spread in a large number of versions over the atomic collisions in solids. It is applied for the first time to the field of ion lithography, in (Mladenov et al., 1985).

We develop a program version, named TRIM-MV (Vutova & Mladenov, 1992), based on the famous Biersack model, to model the ion scattering in multilevel multicomponent amorphous targets. Depending on the different penetration velocities an accelerated ion changes its own charge state being neutralized at the sample surface or going to the highest charged state during its movement through the irradiated sample. The critical velocity for that change is the average orbital velocity of the electron in the statistical atom model of Thomas-Fermi:

where Z_{1} is the atomic number of the penetrating particles. Ions penetrating in the solid with velocity less than V_{cr}, move as neutral atoms, taking an electron from the target material. In the case of V>V_{cr}, they loss own electrons. At these higher energies the interactions (and energy losses) due to collisions with the target electrons predominate while at the lower energies the interactions with the atom nucleus lead to more important type of energy losses. In the common case the total ion energy loss is calculated as a sum:

The nuclear and electronic stopping powers are assumed to be independent of each other. The energy losses for different energies and combinations of types of penetrating ions and target atoms are known. Usually, these energy losses are presented as universal functions of dimensionless parameters ρ and ε (for the ranges and the energy, respectively). These parameters are connected with the real ion range R and the energy E by the relations:

Here M_{1} and M_{2} are the atomic masses, Z_{1} and Z_{2} are the atomic numbers of the bombarding particle and the target atom; N is the density of scattering atoms, a is the shielding parameter of the interaction atomic potential. The Molier potential is appropriate approximation at low energies of the penetrating particle:

where the shielding parameter *a* is presented by the Firsov distance:

The ranges through the penetrated particles’ trajectories are defined by the total ion energy losses and are calculated using the relation:

In the case of low-energy the electronic energy losses are calculated on the basis of the Lindhard and Scharff formalism (Biersack & Haggmark, 1980):

and *к* is given by:

where *ζ*
_{
е
} is a constant the value of which is about *Z*
^{
1/6
}.

The electronic energy losses for high energies are calculated using the Bethe-Bloch equation as Z^{2} is replaced by Z_{1}Z_{2}. The following interpolation expression is also used:

The nuclear losses are evaluated using the transferred (to the target atom) energy T at an elementary collision:

For calculation of the scattering angles θ at the nuclear collisions, using a “center of mass” coordinate system, the impact parameter value p is determined (the laboratory coordinate system and the “center of mass” system are identical in the case of collisions of the penetrating electrons with the target atoms, while these systems are different for penetration of ions in solids). The impact parameter is determined using random numbers:

Then the value of the scattering angle is calculated using the Biersack and Haggmark “scattering triangle” (Fig.19):

where *R*
_{
1
} and *R*
_{
2
} are the trajectory curvature radiuses measured in the point of the minimum distance between the colliding particles *M*
_{
1
} and *M*
_{
2
}; *r*
_{
0
} is the distance between the particles at the moment of minimum distance between them. The corrections *δ*
_{
1
} and *δ*
_{
2
} are low and have smaller values. The values of *r*
_{
0
} are determined using the value of the penetrating ion energy (in the center of mass system):

The values of R_{1} and R_{2} are calculated using (Biersack & Haggmark, 1980):

where:

The values of *δ*
_{
1
} and *δ*
_{
2
} depend on the values of the ion energy E and of the target atomic mass (i.e. r_{0}). In (Biersack & Haggmark, 1980) their values are given in the case of the Coulomb potential using an appropriate analytical expression. For the Moliere potential the values of *δ*
_{
1
} and *δ*
_{
2
} are evaluated using the Robinson values of the scattering integral. The values of these corrections come to values of the scattering angle identical with the Rutherford scattering angle.

After evaluating the value of the scattering angle θ in the “center of mass” system, we use the following formula to calculate the value of θ_{ls} in the laboratory system:

In order to achieve high computer efficiency, a suitable expanding of the distances between consecutive collisions at high energies is applied. The following expression is used instead of Eq.(41):

where the function *L(Е,М*
_{
1
}
*,М*
_{
2
}
*)* is evaluated as in (Wilson et al., 1977).

As in the case of simulation of electron scattering in the target, the penetration of each particle from one statistically significant extract (about 1000 – 10 000 particles) is considered as an elementary collision sequence with the target atoms. Characteristic changes in the particle motion for each collision, assuming a straight line trajectory between two collisions are calculated. The calculation continues until the cut-off conditions are not reached (the value of the penetrating particle is less than 100 eV or the particle is out of the sample). Bragg’s rule (Biersack & Haggmark, 1980) is assumed to be valid. The role of the displaced secondary atoms in cascades is taken into account. The program involves the same MC technique to calculate the kind of the atom taking part in the collision (Eqs.(5-6)). The value of the free path between the collisions is also given by Eqs.(7-8). The difficulties concerning calculation of values of the different angles and the corresponding projection of the ranges (the straight line free paths between the sequential collisions of the projectile ion with two of the resist atoms) are overcome in the same way using Eqs.(11-18). The low-energy nuclear reactions (or anomalies) are not taken into account due to their minor influence on the ion stopping power.

### 7.2. Results and discussion

It is necessary to understand deeply the mechanisms of the processes at IBL in the critical dimension region below 150 nm and to predict by simulation more exactly the ion penetration, the absorbed energy dispersion and the resist modification as well as the subsequent process of the image development. For this aim more details of the calculated results can be experimentally confirmed or rejected. Accumulation and analysis of the calculation results after various simulation stages and of the physical parameters included in the mathematical modeling in nano-image region of IBL is also a question of present interest. The computer simulation allows us to obtain more numerical data and to reveal the dependence of these values on the incoming ion number or mass, as well as on ion energies.

The values of projected and lateral ranges as well as the mean electronic and nuclear stopping powers have been calculated for various penetrating particles at different energies (Mladenov et al., 1985; Vutova & Mladenov, 1992; Vutova & Mladenov, 1994; (a) Vutova & Mladenov, 2001; Mladenov et al., 2001; Vutova, 2007), using our simulation tool TRIM-MV.

Fig.20 presents simulated trajectories for He ions in PMMA. The dependence of the electronic energy losses on the penetration depth in the PMMA resist is shown in Fig.21. It is evident that the contribution of the recoil particles to the total electronic losses is negligible. The nuclear energy losses are presented in Fig.22. We distinguish nuclear collisions accompanied with an energy transfer exceeding 25 eV (curve 1) and collisions with an energy transfer less than 25 eV (curve 2). The former is just the fraction of the nuclear energy losses which can be stored as defects in the material while the latter is the energy which can be converted directly into random (thermal) nuclear motion. The division line in the transferred energies, here chosen as 25 eV, for the creation of permanent defects on one side, and thermal vibrations on the other, is not well known for such polymers as PMMA. According to the Lindhard, Scharff and Schiott (LSS) theory, a comparison between Fig.21 and Fig.22 confirms that the electronic stopping is the dominant form of energy losses during most of the particle path through the solid, for the energies and light masses considered here. Only when the particle is about to come to rest, the nuclear energy loss will play a significant role. From Fig.22 one can also conclude that the nuclear energy losses, for which the transfer energies exceed the threshold value (here 25 eV), are more important than the energy losses due to the collisions with transfer energies less than 25 eV.

Depth distributions of the electronic and nuclear energy losses are presented in Fig.23. The lateral spread is low for the first few thousands Å in depth and then increases to a maximum value of about two thousands in the last part of the ion range – Fig.23. It can be seen that the lateral spread of the nuclear stopping is more pronounced than the electronic stopping. Figure 24 presents the two-dimensional range distribution of 60 keV implanted ^{4}He ions.

Fig.25 shows the projected penetration ranges R_{p}= < t_{x}> (the mean value of the projections t_{x} of the ranges t on the axis x coincident with beam axis and with the sample depth; the brackets < > mark the average procedure) and the longitudinal (penetration) range deviation ΔR_{p}= <(t_{x} - <t_{x}>)^{2}>^{1/2} of the ions ( H^{+}, He^{+}, B^{+}, Ar^{+} and Ga^{+}) in a wide energy region.

In Table 3 are shown the calculated projected penetration ranges *R*
_{
p
}= *<t*
_{
x
}
*>* and *ΔR*
_{
p
}
*= <(t*
_{
x
}
*- <t*
_{
x
}
*>)*
^{
2
}
*>*
^{
1/2
} of the ions ( H^{+}, He^{+}, Li^{+,} Be^{+}, B^{+}, Ar^{+} and Ga^{+}) with an energy of 100 keV in PMMA. Additionally, Table 3 shows the transverse ion range deviation *ΔY = <(t*
_{
y
}
*- <t*
_{
y
}
*>)*
^{
2
}
*>*
^{
1/2
}
*= <(t*
_{
y
}
*)*
^{
2
}
*>*
^{
1/2
} (the projected length *t* of the penetrated ion ranges between two collisions on the direction axis *y* parallel to the sample surface, characterizing the width of the irradiated line) in the case of mentioned irradiating ions.

From these data it is concluded that only H^{+}, He^{+} and B^{+} cover the requirements of the present day IBL which aims to achieve 150 nm resolution patterns because the values of *R*
_{
p
} are sufficient high at chosen resist thickness and these ions have appropriate values of *ΔY* for specified resolution. Among them He^{+} ions are chemically inactive and lighter than the other ions, except than H^{+}, and therefore present some advantages. The ranges of heavier ions Ar^{+} and Ga^{+} at the mentioned accelerating voltages are satisfactory only for the resist thickness up to 100-150 nm. Note that heavier ions (Ar^{+} and Ga^{+}) can be used for higher resolution. The issues are Coulomb interaction and the mentioned limited penetration depth. The latter issue can be solved with multilayer resist with top-surface imaging.

Ion | Mi | Zi | Rp [nm] | ΔRp [nm] | ΔY [nm] | (dE/dx)el [eV/nm] | (dE/dx)n [eV/nm] |

H+ | 2 | 1 | 994 | 58 | 107 | 120 | 0.13 |

He+ | 4 | 2 | 796 | 78 | 115 | 185 | 1.96 |

Li+ | 6 | 3 | 720 | 88 | 120 | 190 | 4.53 |

Be+ | 9 | 4 | 562 | 80 | 108 | 195 | 10.16 |

B+ | 11 | 5 | 514 | 77 | 91 | 205 | 16.64 |

Ar+ | 40 | 18 | 183 | 47 | 30.9 | 371 | 158 |

Ga+ | 70 | 31 | 140 | 35 | 22 |

One can see (Table 3) that the values of the transverse range deviation *ΔY* of the studied light ions are between 90 nm and 120 nm. At these values of the ion ranges *ΔY*, the patterning is affected by the ion scattering. The changes of the latent image dimensions in comparison with those of the irradiated image are considerable for high resolution layouts with critical dimensions of 50-150 nm. The situation is the same as in the case of electron lithography where the proximity effect and the related corrections of exposed patterns take part. In Table 3 calculated data of the energy losses for these light ions in PMMA are shown also where: *<dE/dx>*
_{
el
} stands for the mean electronic energy losses and *<dE/dx>*
_{
n
} stands for the mean nuclear energy losses. Data calculated for the all other moments characterising the spatial distribution of the He^{+} with energies in the range of 30-300 keV penetrating in PMMA are shown in Fig.26-27. The upper curve (Fig.26) represents the calculated ranges through the penetrated particle trajectories R_{l} = <t>. The mean quadratic deviation ΔR_{l} = <(t - <t>)^{2}>^{1/2} and ΔY are also shown in Fig.26.

The curve situated under the horizontal axis in Fig.27 represents the moment of the third order characterising the distribution symmetry i.e. the skewness γ=<t^{3}
_{x}/ΔR_{p}
^{3}>. For example: γ=0 means that the distribution of the ranges is symmetrical. The next moment of the projected range distribution is the kurtosis β=<t_{x}
^{4}/ΔR_{p}
^{4}>, which characterises the tail of the projected range distribution (curve 3 in Fig.27). For example: β=3 corresponds to the normal distribution; β=6 is in the case of an exponential distribution; the higher values of β mean that longer distribution tails exist. Curve 1 shows the transverse kurtosis concerning the distribution symmetry in a width. Curves 2 and 4 represent the moments <t_{x}t_{y}
^{2}/t_{y}
^{3}> and <t_{x}
^{2}t_{y}
^{2}/t_{y}
^{4}>, respectively. By such calculations one can see that the scattering and lateral dispersion of the penetrating particles through a resist are small for the ions, which are three orders of magnitude heavier than the electrons. Therefore the proximity effects that are the major problem in the EBL are minimal in the IBL.

Fig.28 shows calculated energy losses in the case of He^{+} bombardment of PMMA where <dE/dx> is the mean electronic energy loss. Indexes 1 and 2 correspond to the primary and secondary ions, respectively. The upper curve represents the electronic energy losses of primary ions, whereas the lower curve describes these losses for secondary ions. The total electronic energy losses (the sum of these losses of primary and secondary ions) coincides with the upper curve (Fig.28) since secondary ions have a small contribution to the total electronic energy losses. The index f of the mean nuclear energy losses <dE/dx>_{nf} means losses of energy through the nuclear collisions (Fig.29). The maximum value of these losses is assumed to be equal to 20 eV. It is assumed that these collisions create only phonons (i.e. heating). The index d of <dE/dx>_{nd} indicates the higher portions of lost energy (>20 eV) through the nuclear type of collisions, which are able to change the resist structure completely by the creation of defects – curve <dE/dx>_{nd,}
Fig.29. The total nuclear energy losses <dE/dx>_{nΣ} are also shown in Fig.29. A comparision between the electronic and nuclear stopping power values in the case of He^{+} in PMMA confirms the more general conclusion, based on the experimental data, that the electronic stopping losses determine the solubility changes in the resists (Mladenov & Emmoth, 1981). As an additional reason can be indicated the fact, that the values of *<dE/dx>*
_{
el
} are more than two orders of magnitude higher than *<dE/dx>*
_{
n
} at the studied kinds of ions and their energies.

Using these data it can be concluded, that the radiation changes due to atom mixing of resist atoms at nuclear collisions of He projectiles with resist components in the studied energy range have a low probability. The situation is opposite to the case of 120 keV Ar^{+} irradiation ((b) Vutova & Mladenov, 2001) where lower solubility rate at higher doses of ion exposure had been explained by ion irradiation defects. In the same time the electronic energy losses are not so high, that erosion, decomposition and change of physical properties of the irradiated resist take place (Braun et al., 1983; Emmoth & Mladenov, 1983).

To obtain the absorbed energy space distribution when exposing an arbitrary pattern, using an arbitrary exposure dose distribution, one must take into account the influence of a large number of exposed points. Due to the large number of calculation, the simplified procedure (described above for EBL) is used to calculate the integral space distribution of the absorbed energy. To calculate the δ_{f}, δ_{b} and η_{E} parameter values (Eq.19), the linear approximation along the resist depth is used. In some ion lithography cases (heavy-ion bombardment, high-energy ion beams, etc) it is necessary to use more than two layers along the resist depth (Vutova & Mladenov, 1994). The values of δ_{f}, δ_{b} and η_{E} change linearly among them.

## 8. Conclusions

In this chapter we describe the basic steps when modeling the ion and electron exposure processes in multilevel multicomponent amorphous targets. The modeling of the exposure process in EIL involves the simulation of particle penetration into materials as well as the calculation of the absorbed energy distribution within the targets. The models presented for the exposure process simulation are realized by means of a software package. Using our simulation tools the effects of: irradiated particles (electrons and ions), exposure doses, substrates and dense patterns are discussed in order to extract the necessary values for high resolution patterning.

EBL simulation of the exposure process was carried out in the case of semiconducting and superconducting substrates using the developed models and the computer programs. For the calculation of energy deposition in the resist film due to a point beam (EDF(r,z)) a MC simulation was applied. The problem concerning the insufficient statistics of the discrete data for the absorbed energy in the case of large lateral distances is overcome. In the Si case, the EDF(r) is approximated by a sum of two Gaussian functions, while in the case of YBCO layers deposited on STO or MgO substrates, the EDF(r) could be approximated as a combination of a double Gaussian and an exponential function. In this case the division of scattering electrons into forward and backward is not the same as for Si samples since the backscattering from YBCO film is important in the middle region of the EDF(r) (from ~ 40 nm to 2 μm). Data for the scattering parameters, such as standard deviations of EDF(r), was obtained. A methodology and a computer program are proposed for the calculation of the absorbed energy integral space distribution when exposing real topological structures in EBL, using analytical expressions. The simulation results were compared to the experimental and other calculated data in the case of Si and YBCO/STO(MgO) substrates and the agreement was good.

A program version, named TRIM-MV, for exposure process modeling at IBL, based on the famous Biersack model, is developed. From the detailed simulation study, significant information on the energy losses, ranges and moments of the space distributions for different ions is obtained. It is concluded that the exposure with He ions with energies from 60 to 100 keV is favorable for sub-150 nm region. On the basis of our experimental and calculated results, it is assumed that the electronic energy losses are responsible for solubility modification of polymer resists in IBL. The nuclear stopping power has an effect on the solubility modification in the case of heavy low-energy ion irradiation at non-linear (multiciphered) development rates.

The obtained results and comparisons show that the models for exposure process modeling at EBL and IBL are adequate and have a very good potential for application in the case of multiplayer sub-quarter-micron patterns. The simulation results are very useful when optimizing particular technological processes for complex layouts in the sub-quarter-micron lithography. The ion lithography should be considered as one of the competitors for further improvement of the integration level of the produced integral circuits.