Power law nature in electron solid interaction

Monte carlo simulation of paths of a large number of impinging electrons in a multi-layered solid allows to define area of spreading electrons (A) to capture overall behavior of the solid. This parameter 'A' follows power law with electron energy. Further, change in critical energies, which are minimum energies lost corresponding to various electrons, as a function of variation in lateral distance also follows power law nature. This power law behavior could be an indicator of how strong self-organization a solid has which may be used in monitoring efficiency of device fabrication.


I. INTRODUCTION
The interaction of beam of energetic electrons with the target solid material technique, which is electron beam lithography, is of great interest in probing material properties (chemical, electrical, physical etc) [1] at sub-micron and nanoscale level [2], and has many applications, modeling radio-induced cellular damages [3], in surface science technology [4], nanolithography techniques [5], fabrication of fractal surfaces [6], various biological applications [7] etc. The impinging energetic electron suffers random collisions from a number of scattering centres with random distribution of potentials in the centres in the solid materials of various layers [3], and the electron follows a stochastic path inside the solid material [1]. The analysis of the electron paths could highlight some of the important properties of the material which will be used in various device fabrications and various other applications.
One of the most important properties of real networks, ranging from social to biological protein-protein networks [8], is power law nature of the network distribution [9] which could be a reflection of fractal nature of the system [10]. Since fractal behavior of the system can be used as an indicator of self-organization in the system [11], one can use this property to identify important patterns and their origin in the system [12][13][14]. The path of the penetrating electrons in solid system is the reflection of organization of the regular scattering centres with random potential distributions, one can identify probable parameters to capture patterns of organization in the solid material. In this work, we try to search for possible parameters to characterize fractal nature of the solid systems using Monte carlo simulation procedure which could be used for various fabrication techniques. In the section II the detailed Monte carlo procedure of electron-solid interaction is described. Simulation results are described in section III, and some conclusions are drawn based on the simulation results. * Electronic address: brojen@jnu.ac.in (Corresponding author)

II. ELECTRON SOLID INTERACTION MODEL
The path traversed by impinging energetic electrons in solid is based on the electron transport within the stochastic formalism of scattering process of electrons with the solid along their trajectory [3,15]. The penetrating electrons encounter randomly distributed scattering centres within the electron interaction range [3], and the electrons undergo complicated brownian paths inside the solid [16]. These electrons with energy (E) move in straight lines between any two scattering centres, and once they suffer interaction with scattering centres of the solid, the change in their directions are defined by (E, θ, φ), where, θ and φ are scattering and azimuthal angles respectively. The solids could be single or multilayered thin film with different distributions of scattering centres in different solid layers.
We consider impinging electrons suffer elastic collisions from scattering centres distributed in the single or multilayered thin solid film, where, differential cross section can be described by classical screened Rutherford's formula [16], where, e is electronic charge, Z is the atomic number of the material and dΩ is solid angle. β is screening parameter to account for electrostatic screening of the nucleus by the orbital electrons, and is given by Thomas-Fermi model of the atomic field [17], where, a o is Bohr radius, is Planck's constant and p is the electron's momentum. The scattering angle θ can be calculated by evaluating total elastic cross section σ = dΩ from equation (1), where, F (θ) is accumulated function of scattering probability [16] which is a function of theta only. In the Monte carlo simulation procedure, θ can be obtained by generating a uniform random number R 1 in [0, 1], and azimuthal angle φ by generating another uniform random number R 2 and using, Thus, θ and φ can be calculated by generating two sets of uniform random numbers, and using equation (3) and (4).

A. Energy loss calculation of traversing electron
The energy of the electron, suffering interaction from the scattering centres distributed randomly along its path, continuously looses its kinetic energy and can be calculated using Behte's continuous slowing down approximation model [19]. This model is a good emperical method for high energetic electrons as compared to ionization energy J i.e. for E >> J, but suffers problem for E ≤ J [20]. However, the model was generalized for all range of energies [21], where, the energy loss ∆E of the penetrating electron a path length L along its trajectory is given by, where, M is the atomic weight of the target material, ρ is the density, N A is the Avogadro's number and C is a constant (C → 1; C < 1). The mean ionization energy J can be obtained from the emperical formula [22], where, J Z → 9.76 as Z → ∞. Further, the sensitivity of the J in Monte carlo simulation can be controlled by taking logarithm of this parameter.

B. Modeling electron path in multi-layered system
The mean free path for single layered system, calculated using equation (1), can be extended for multilayered system by defining a probability P m (u) that the electron once scattered from first layer is not scattered until mth layer [23,24], and can be obtained from the following equation, where, Γ m is scattering probability per unit length in mth layer. Boundary condition is taken as P 1 (0) = 1. Solving equation (7) one can arrive at P m+1 (u − u m ) = P m (u), where, u m is the distance between mth and (m + 1)th layers along z-axis. This P m (u) can be related to F (u) by, F (u) = 1 − P m (u) = 1 − R 1 , and θ can be solved using equation (3). The mean free path for single layer system is calculated as λ 1 = ∞ 0 uP 1 (u)du = −1/Γ 1 , where, P 1 (u) = exp(−uΓ 1 ). Similarly, mean free path for two layered system is given by, . Proceeding in the same way, mean free path for m layered system can be calculated using, From equation (8) one can able to calculate u m of impinging electron in mth layered material system. Now, starting from an initial vector (x 0 , y 0 , z 0 ) T , we can trace the path of the penetrating electron in m-layered system using the following recursive procedure, Thus the path of the penetrating electron in m-layered material system can be traced using the recursive procedure in (9) with average energy loss ∆E within the Monte carlo simulation.

III. RESULTS AND DISCUSSION
We consider single layered GaP and double layered system of GaP with resist P M M A (C 5 H 8 O 2 ), and simulated using the Monte carlo simulation procedure described in the previous section to trace the trajectory of impinging energetic electrons in the systems and energy loss. The initial positions of the penetrating 500 electrons are taken as the same as (x 0 , y 0 , z 0 ) T = (0, 0, 0) T for a range of energy keV. The trajectory of each electron of energy E is calculated using Monte carlo procedure (9), and energy loss (∆E) during the process of penetration is obtained using equation (5). In our simulation the ionization energies of C, H and O are taken to be 78eV, 18.7eV and 89eV, and for GaP , equation (6) is used to calculate its J. Since the ionization energy of both GaP and P M M A are of the order of eV, E >> J in our case and therefore, we take C < 1 in equation (5) during the calculation. The thickness of the P M M A resist is taken to be 10 −6 metre in the double layered calculation.
A. Power law nature of electron spreading area The paths of the impinging electrons (500 electrons' paths) in single layer GaP system with different energies ([2-100]keV) are calculated ( Fig. 1 left panels) and two dimensional areas of the spreading electrons (circles in the figures) for different energies are obtained. The area of each circle is calculated for ten ensembles, and average of minimum and maximum areas bounded to the two dimensional electron spreading area of 500 electron trajectories (error bars in the plot in middle lower panel of Fig. 1). This calculated areas A as a function of E is found to follow the following power law behaviour, The straight line is the fitting curve to the calculated data. The value of γ is found to be 1.98. We then calculated the As in double layered material system (1µm thickness for P M M A and rest for GaP ) for various electron energies [2-100]keV ( Fig. 2 left panels). The behavior of A with respect to E (Fig. 2 middle upper panel) has three regimes, left (electrons paths are within single first layer only) and right (dominated by second layer as compared to first layer) regimes follow similar power law given by equation (10), and the values of γ are found to be 2.1 and 1.97 respectively. The middle regime, which is due to contributions from both first and second layers, and does not follow exact power law nature.

B. Power law behaviour in energy loss
The energy loss of impinging 1500 electrons in single layer system are calculated as a function of lateral distance R of the electron trajectories for different energies ( Fig. 1 right panels). The values of E sharply drop after a certain value of R for different values of E showing that the electrons do not have sufficient energies to penetrate further in the solid. We then calculated these critical E c and R c for different Es in the range [10-100]keV. Then we calculated possible changes in these critical energies (∆E) as a function of ∆R starting from lowest (E c , R c ), and error bars are standard deviations of the thicknesses of the drop curves (lines drawn parallel to E axis in the Fig. 1 right panels). The calculated ∆E again follows power law with ∆R as follows,