Three-Dimensional Simulations with Fields and Particles in Software and Inflector Designs

Particles and fields represent two major modeling paradigms in pure and applied science at all. In this paper, a methodology and some of the results for three-dimensional (3D) simulations that include both field and particle abstractions are presented. Electromagnetic field calculations used here are based on the discrete differential form representation of the finite elements method, while the Monte Carlo method makes foundation of the particle part of the simulations. The first example is the simulation of the feature profile evolution during SiO2 etching enhanced by Ar + /CF4 non-equilibrium plasma based on the sparse field method for solving level set equations. Second example is devoted to the design of a spiral inflector which is one of the key devices of the axial injection system of the VINCY Cyclotron.


Introduction
As already pointed out, fields and particles are two the most important modeling paradigms in pure and applied science at all.Particles are primarily used in one of the two ways in large scientific applications.The first one is to track sample particles to gather statistics that describe the conditions of a complex physical system.Particles of this kind are often referred to as "tracers".The second one is to perform direct numerical simulation of systems that contain discrete point-like entities such as ions or molecules.In both scenarios, the application contains one or more sets of particles.Each set has some data associated with it that describes its members' characteristics, such as charge, mass or momentum.The state of the physical system is defined by these data [1].
Particles typically exist in a spatial domain, and they may interact directly with one another or with field quantities defined on that domain.A field, on the other hand, defines a set of values on a region of space.In order to specify a field one must specify the locations at which the field's values are defined, and describe what happens at the boundaries of that region in space.For this purpose we use meshes which are discrete representations of the physical domains.Particle abstractions are usually designed to be used in conjunction with fields.Some types of interpolators are used as the glue that bind these together, by specifying how to calculate field values at particle (or other) locations that not happen to lie exactly on mesh points.Interpolators are used to gather values to specific positions in a field's spatial domain from nearby field elements, or to scatter values from such positions into the field.
In this paper we present two examples of the simulations of this type, different in their origin but sharing the same software design concept and implementation components.The first is from the area of the applications of the plasma technologies in microelectronic circuits manufacturing.The second is a part of the ongoing work on the design of the spiral inflector for the VINCY Cyclotron.

Software Design Issues
High-performance scientific applications are notoriously difficult and expensive to develop and maintain.Very little progress had been made in scientific programming productivity over decades.On the other hand, in the business programming community giant steps have been made in productivity, through software technology including component-based programming and frameworks.Attempts to transfer this technology to the scientific community have largely been unsuccessful until recently.However, recent advances in software technology signify that the hard to get goal of increased high-performance scientific programming productivity is possible.
Although Object-Oriented Programming itself was not the main reason for the improved productivity the business community (the rise in productivity comes from the use of Object-Oriented Design and Component-Based Programming), the choice of the programming language is of the utmost importance.Development of complex scientific applications requires GUI toolkits, numeric libraries, graphical facilities, libraries for interfacing with code written in other languages, etc.Many scientific programmers and researchers still dream of a general-purpose language that is so expressive, elegant, and efficient that it can cover all needs without significant inconveniences.However, no current language can satisfy all these requirements natively, and there is no a research language that could possibly do that in the short-to-medium term.
Our simulation framework (written in C ++ ) is based on Fltk library [2] for building GUI elements, Vtk toolkit [3] for 3D visualization needs, and Pthreads library [4] for managing multiple execution threads.Additional libraries and toolkits are used depending on the specific requirements of the particular applications.

Electromagnetic Field Calculations
In this approach the scalar electrostatic potential is a 0form, the electric and magnetic fields are 1-forms, the electric and magnetic fluxes are 2-forms, and the scalar charge density is a 3-form.The basic operators are the exterior (or wedge) product, the exterior derivative, and the Hodge star.Precise rules (i.e. a calculus) prescribe how these forms and operators can be combined.In this modern geometrical approach to electromagnetics the fundamental conservation laws are not obscured by the details of coordinate system dependent notation.By working within the discrete differential forms framework, we are gauranteed that resulting spatial discretization schemes are fully mimetic [5,6].
The calculations of the electric fields in this paper are performed by integrating a general finite element solver GetDP [7,8] in our simulation framework.GetDP is a thorough implementation of discrete differential forms calculus, and uses mixed finite elements to discretize de Rham-type complexes in one, two and three dimensions.Meshing of the computational domain is carried out by TetGen tetrahedral mesh generator [9].TetGen generates the boundary constrained high quality (Delaunay) meshes, suitable for numerical simulation using finite element and finite volume methods.As a part of the post-processing procedure the electric fields, obtained on the unstruc-tured meshes, are recalculated on the Cartesian rectangular domains containing the regions of the particles movement.In this manner, the electric field on the particles could be calculated by simple trilinear interpolation.

Simulation of Etching Profile Evolution
Refined control of etched profile in microelectronic devices during plasma etching process is one of the most important tasks of front-end and back-end microelectronic devices manufacturing technologies.The profile surface evolution in plasma etching, deposition and lithography development is a significant challenge for numerical methods for interface tracking itself.Level set methods for evolving interfaces [10,11] are specially designed for profiles which can develop sharp corners, change topology and undergo orders of magnitude changes in speed.They are based on a Hamilton-Jacobi type equation [12] for a level set function using techniques developed for solving hyperbolic partial differential equations.A simple model of Ar + /F etching process that includes only pure chemical and ion-enhanced chemical etching mechanisms is presented in details.

Sparse Field Level Set Method
The basic idea behind the level set method is to represent the surface in question at a certain time t as the zero level set (with respect to the space variables) of a certain function φ(t, x), the so called level set function.The initial surface is given by {x|φ(0, x) = 0}.The evolution of the surface in time is caused by "forces" or fluxes of particles reaching the surface in the case of the etching process.The velocity of the point on the surface normal to the surface will be denoted by V(t, x), and is called velocity function.For the points on the surface this function is determined by physical models of the ongoing processes; in the case of etching by the fluxes of incident particles and subsequent surface reactions.The velocity function generally depends on the time and space variables and we assume that it is defined on the whole simulation domain.At a later time t > 0, the surface is as well the zero level set of the function φ(t, x), namely it can be defined as a set of points {x∈ℜ n|φ(t, x) = 0}.This leads to the level set equation in the unknown function φ(t, x), where φ(0, x) = 0 determines the initial surface.Having solved this equation the zero level set of the solution is the sought surface at all later times.Actually, this equation relates the time change to the gradient via the velocity function.In the numerical implementation the level set function is represented by its values on grid nodes, and the current sur-face must be extracted from this grid.In order to apply the level set method a suitable initial function φ(0, x) has to be defined first.The natural choice for the initialization is the signed distance function of a point from the given surface.This function is the common distance function multiplied by −1 or +1 depending on which side of the surface the point lies on.As already stated, the values of the velocity function are determined by the physical models.In the actual numerical implementation Equation ( 1) is represented by the upwind finite difference schemes (see ref. [10] for the details) that requires the values of this function at the all grid points considered.In reality the physical models determine the velocity function only at the zero level set, so it must be extrapolated suitably at grid points not adjacent to the zero level set.
Several approaches for solving level set equations exist which increase accuracy while decreasing computational effort.They are all based on using some sort of adaptive schemes.The most important are narrow band level set method [10,11], widely used in etching process modeling tools, and recently developed sparse-filed method [13], implemented in medical image processing ITK library [14].The sparse-field method use an approximation to the distance function that makes it feasible to recompute the neighborhood of the zero level set at each time step.
In that way, it takes the narrow band strategy to the extreme.It computes updates on a band of grid points that is only one point wide.The width of the neighborhood is such that derivatives for the next time step can be calculated.
Here we will present some calculations illustrating our approach for etching profile evolution simulation.All the calculations are performed on 128 × 128 × 384 rectangular grid [15].The initial profile surface is a rectangle deep with dimensions 0.1 × 0.1 × 0.1 μm.Above the profile surface is the trench region.From its top the particles involved in etching process come from, while bellow is the non-etched material.The actual shape of the initial surface can be described using simple geometrical abstractions.In the beginning of the calculations this description is transformed to the initial level set function using fast marching method [10].The evolution of the etching profile surface with time is shown in the following figures.In Figure 1(a) the results obtained for a test calculation performed with constant velocity function V = V 0 = 5 nm/s (purely isotropic etching case) are shown.It is supposed that only the bottom surface could be etched; i.e. that the top and the vertical surfaces belong to photo-resist layer.Behavior of the etching profile is as expected.
In the case of the ion enhanced chemical etching the dependence of the surface velocity on the incident angle is simple [16]: V = V 0 cosθ.The pure chemical etching velocity, or more precisely the etching yield, does not depend on the incident angles.This case can be safely treated by the upwind scheme and using the Lax-Friedrichs scheme would lead to unnecessary rounding of the profile surface.The high aspect ratio (depth/width) etching is a common situation in semiconductor technologies.The evolution of the etching profile, when etching rate is proportional to cos(θ), is presented in Figure 1(b).This is the simplest form of angular dependence, but it describes the ion enhanced chemical etching process correctly [16].In this case we expect that the horizontal surfaces move downward, while the vertical ones stay still.This figure shows that it with optimal amount of smoothing gives minimal rounding of sharp corners, while preserving the numerical stability of the calculations.Actually, this is one of the most delicate problems in the etching profile simulations.
A comprehensive simulation of etching requires knowledge of the etching rates at all the points of the profile surface during the etching process.These rates are directly related to fluxes of the etching species on the profile surface, which are themselves determined by the plasma parameters in the etching device.Electrons do not contribute directly to the material removal, but they are the source, together with positive ions, of the profile charging that has many negative consequences on the final outcome of the process especially in the case of insulating material etching, SiO  fluxes incident on the profile surface.Our simulation concept [20] is similar in spirit to the 2D simulations presented in [21,22] and especially [23], where charging effects in 3D rectangular trench were analyzed.Monte Carlo technique is the only feasible method for calculating particle fluxes in 3D geometries.Trench wall charging strongly influences the charged particles motion and, consequently, particle fluxes which themselves determine the local etching rates.
Since the trench boundaries have no regular (rectangular) shape in our simulation, finite element calculations was used for the calculation of the electric field.As the etching profile is not known in advance (it is a result of the calculations itself), the problem of meshing is extremely difficult.In Figure 3 the electrostatic potential map is shown for a test case calculation, for illustrating purposes only.Electric field obtained in that way is used in standard leap-frog particle moving scheme.

Inflector Design
This part work is related to the activities concerning the design of the spiral inflectors for this facility, which is the second example of coupled field and particle simulations treated in this paper [24].The main roles of the inflector are to bend the injeced ion beam from its initial direction onto the cyclotron's madian plane and, together with the optical elements in the transport line, to match the beam emittance to the cyclotron's acceptance.The electrostatic inflector consists of two biased and two grounded electrodes.The electric field produced by these electrodes exerts a force on the ions, simultaneously bending and focusing the beam.The design of an electrostatic inflector is complicated by the fact that in addition to the electrostatic force produced by its electrodes, the ions are also subjected to a magnetic force produced by the simultaneously bending and focusing the beam.The design of an electrostatic inflector is complicated by the fact that in addition to the electrostatic force produced by its electrodes, the ions are also subjected to a magnetic force produced by the magnetic field near the center of the cyclotron.This effect must be taken into account in designing the inflector.Different types of inflectors have been devised [24][25][26][27] for inflecting the axially injected ion beam into the cyclotron median plane.In modern, variable energy, multi-particle, compact cyclotrons, the minimal gap between the magnetic poles tends to be very small (few centimeters-to provide high flutter and high magnetic circuit efficiency).This fact imposes severe restrictions on inflector dimensions, as well as specific demands concerning its optical properties.Owing to its flexibility and relatively low voltage needed for its operation, the electrostatic spiral inflector has become widely used in the multi-particle compact cyclotrons.Due to the complexity of the shape of its electrodes, this particular type of inflector has been chosen as illustrative (see Figure 4).
The ion trajectories are calculated by solving Newton-Lorenz equations of motion using Boris integrator scheme [1], since magnetic field must be taken into account also.
In Figure 5 an electrostatic potential map at the exit plane of the inflector is shown.The voltages applied to the upper and the lower electrodes are taken to be −1 V and 1 V in the potential calculations.

Conclusion
In this paper we presented some results of two different types of combined particle-field simulations (based on the unified software framework), in which the authors are involved.The first is connected with the applications of the plasma etching technologies in microelectronic circuits manufacturing.Obtained results show that more realistic calculations should include the effects of polymer deposition, profile charging, more complex surface reactions set, as well as better statistic (greater number of particles) in the Monte Carlo step of the calculations, which is now limited by the available computational resources.Another is a part of the ongoing work on the design of the spiral inflector for the VINCY Cyclotron.In present phase this simulations are preliminary, and more comprehensive analysis that includes the detailed maps of central region magnetic field, particle bunches and interparticle interaction is necessary.
2 for example.The energy and angular distribution functions for Ar+ ions (IEDF, IADF) and electrons (EEDF, EADF) are shown in Figure 2.They are obtained by particle-in-cell (PIC) calculations using XPDC1 code [17-19].These data are used as the boundary conditions for the calculations of ion Copyright © 2013 SciRes.JSEA Three-Dimensional Simulations with Fields and Particles in Software and Inflector Designs 393

Figure 3 .
Figure 3.An example of the electrostatic potential map of the charged feature profile.

Figure 5 .
Figure 5. Electrostatic potential map at the exit plane of the spiral inflector.