Computer Simulation of Metal Ions Transport to Uneven Substrates during Ionized Plasma Vapour Deposition

We present a computational study of processes taking place in a sheath region formed near a negatively biased uneven substrate during ionized plasma vapour deposition.The sputtered metal atoms are ionized on their way to substrate and they are accelerated in the sheath near the substrate. They are able to penetrate to high-aspect-ratio structures, for example, trenches, which can be, therefore, effectively coated. The main technique used was a two-dimensional particle simulation.The results of our model predict the energy and angular distributions of impinging ions in low-pressure conditions which are characteristic for this method and where typical continuous models fail due to unfulfilled assumptions. Input bulk plasma properties were computed by a “zero dimensional” globalmodel which took into accountmore physical processes important on a scale of the wholemagnetron chamber. Output parameters, such as electrostatic potential, energy of ions, and ion fluxes, were computed for wide range of conditions (electron density and substrate bias) to show the influence of these conditions on observed phenomena, penetration of sheath inside the trench, deceleration of argon and copper ions inside the trench, and local maxima of ion fluxes near the trench opening.


Introduction
When modifying metal or semiconductor solids during various material-treatment technologies, there is an important group of methods connected with both inert gas and chemically active plasmas, oxidation, sputtering, plasma deposition, plasma-enhanced chemical vapour deposition, plasma immersion ion implantation, and so on.In all these cases the detailed understanding of plasma-solid interaction is of a great importance.One of these techniques is the ionized plasma vapour deposition (IPVD) or more specifically Inductively Coupled Plasma Magnetron Sputtering (ICP-MS).
The deposition technique IPVD was presented in [1] as an improvement in filling of high-aspect-ratio features in microelectronics fabrication.Metal atoms sputtered from a target move through inductively coupled plasma of low pressure but high electron density where they are ionized.Ionization degree of metal atoms and flux fraction of metal ions can reach very high values during IPVD.Paper [2] described ionization up to 80% for aluminium atoms.Flux fraction of copper ions exceeding 40% was shown in [3].Afterwards, the metal ions together with ions of the carrier gas move towards the negatively biased substrate.In the vicinity of the substrate, that is, in the presheath and sheath, the ions are accelerated and their angular distribution gets narrower.Ion energy can be regulated by the substrate bias.The flux of metal ions depends on power of the coil which generates the secondary plasma.The flux of ions can be regulated independently on their energy.Accelerated ions can reach the bottom of highaspect-ratio trenches which are hence coated all over their surface [4].Furthermore, the high-energy metal ions and ions of carrier gas are able to resputter the overhanging deposit which is often created at the opening of the trench and which makes further coating impossible [5,6].

Description of the Simulation
In our contribution we studied the physical and chemical processes on surfaces, both flat and uneven, of solids immersed into plasma.The method of solution was the 2 Advances in Materials Science and Engineering computer simulation as an increase in performance of computers in recent years allowed computer modelling to be included among other research tools as a valuable source of information.Plasma is an environment where phenomena over multiple length and time scales are coupled.There exists two basic modelling techniques, the particle models which are able to capture microscopic effects but they demand huge computational resources, while the macroscopic fluid models which are not so computationally demanding however give results with only limited accuracy caused by the lack of microscopic information.To solve specific problems the third way can be used, the hybrid models that take advantage of both fluid and particle modelling approach [7,8].
Particle models, or more precisely PIC-MCC (i.e., Particle-in-Cell, Monte Carlo collisions) models, are quite frequently used techniques to simulate plasma behaviour near immersed solids.For example, similar models were used to simulate conditions during plasma immersion ion implantation process [17,18].Simulation [19] is an example of a common hybrid approach, electrons are treated continuously and ions as particles.Another approach, an iterative hybrid model [20], was used to study problem of plasma diagnostics.We also tried to develop a fully 3-dimensional self-consistent particle model [21] which has no symmetry requirements.However, its computational requirements are very high and hence the size of the mesh and the computational domain are limited.More complex model of IPVD [22] includes processes on the surface, the deposition itself including resputtering and movement of physisorbed particles on the surface.However, the plasma-related part is a continuous model.Our model is intended for precise determination of plasma behaviour near the substrate and it can be further supplemented with the surface-related phase of IPVD using Monte Carlo-based methods.Here, the microscopic nature of particle model results can be of great use.

Obtaining Input Plasma Parameters.
A simple volumeaveraged model was adopted from [23] to estimate the main plasma parameters of argon plasma with addition of sputtered and partially ionized metal, in our case copper.The main processes taken into account are summarized in Table 1 along with the list of their parameters and references from which they have been taken.The rate constants of electron collisions were numerically derived at runtime using cross sections from the literature.The processes mentioned in Table 1 were connected with each other in generationloss balance equations which were solved as a system together with the conservation equation for metal atoms and quasineutrality condition.The equations for individual processes can be obtained from article [10] which deals with aluminium deposition rather than copper.The solution of Ar/Al plasma was also repeated to test the model.

Particle Model.
Dimensionality of the particle model is 23V, that is, two spatial dimensions and three dimensions of velocity.Particles follow two-dimensional trajectories and the electrostatic interaction (Poisson's equation) is solved in 2D.The third dimension of velocity follows Gaussian distribution and plays a part during collisions.
The particle model follows trajectories of three types of charged particles, electrons, argon ions, and copper ions.Their trajectories are influenced by Monte Carlo collisions with neutral argon atoms and an interaction with electrostatic forces originating from the biased substrate and from other charged particles in the computational domain.Copper neutral atoms do not participate in particle model collisions because of their negligible density.In this study we assumed total averaged copper density in the chamber 2 × 10 12 cm −3 which is nearly three orders of magnitude lower than the density of argon.
Instead of random free path, our method of collision simulation generates random time to next collision.The cross section  0 of the imaginary "null collision" satisfies where ] max is the lowest possible constant collision frequency and  is the relative velocity.Unlike the widely used null collision method, this approach does not neglect the velocity of scattering centres.This is important especially for the movement of ions whose velocity is comparable to the velocity of neutral collision centres, at least outside the sheath.This method is derived and described in detail in [24].We included following scattering processes of charged species with argon neutrals: Ar + : elastic scattering and charge transfer [25], electrons: elastic scattering, excitation to the first excited level, and ionization [12], and Cu + : elastic scattering [26].
The electrostatic force was computed in the usual manner using the Particle-in-Cell method.Poisson's equation was solved numerically at each time-step using the sparse system solver UMFPACK [27].For the integration of Newton's equation of motion we used Verlet velocity algorithm [28].The computational domain was square-shaped with side of 5 mm.On one side, there was substrate with the twodimensional trench with size of 2 mm (height) and 0.2 mm (width).The substrate was biased and the given potential, with respect to the plasma potential, was used as Dirichlet boundary condition along the substrate surface.The opposite side of the domain acted as a connection with bulk plasma.Particles could leave the domain through this side and other particles came through it to the domain from a source where we assumed Maxwell velocity distribution of charged particles with corresponding temperatures.The condition  = 0 volts was applied on this side.Remaining sides had periodic boundary conditions.The Cartesian mesh size was set to 1000 × 1000 and on average there were 1 × 10 7 particles in the domain.Figure 1 shows a diagram of the computational domain.

Global Model Results.
For the validation of the volumeaveraged model we used the same input conditions as in [23].Figure 2(a) shows a good agreement of results, electron temperature as a function of total amount of metal in the chamber which contains aluminium.Figure 2(a) also shows an interesting comparison of Al and Cu behaviour in the magnetron chamber.Increasing amount of sputtered metal (caused, e.g., by increasing magnetron power) results in cooling of electrons.This effect is much stronger for Al.Ionization energies are  ,Al = 5.99 eV and  ,Cu = 7.73 eV.There are certainly more electrons in the chamber capable of ionizing Al then Cu and hence they lose ionization energy more easily.It is also in agreement with the deceleration of the temperature decrease.Figure 2(b) shows the development of ionization degree.As long as the electron temperatures are similar for both cases (small amount of metal), aluminium is significantly more ionized.But faster decrease of electron temperature in chamber with aluminium causes also faster decrease of ionization degree.
However, for our simulations we used conditions typical for deposition of copper which were used in the experiment in [3].One of the main variables in this study was RF power.It was shown that the RF power can effectively regulate electron density and consequently the ion-to-neutral flux fraction which is one of the crucial parameters in IPVD process.Therefore, as one of our studied parameters we chose RF power which corresponds linearly to electron density.Other parameters remained the same as in [21].Results of the global model, some of which were used in the particle model, are shown in Table 2.

Particle Model Results. The particle model was run 40x;
for each of these four values of RF power (100 W, 400 W, 700 W, and 1000 W) we simulated 10 different values of substrate bias (−20 V, −40 V, . .., and −200 V).At first, we checked the shape and size of the sheath.Figure 3 shows two contour lines of electrostatic potential, 0.975 ×   and 0.875 ×   .The substrate bias   was −40 V in panel (a) and −200 V in panel (b).In both panels two different RF powers are compared.As expected, the sheath is thicker for lower RF power (electron density) and for higher substrate bias.The penetration of sheath into the high-aspect-ratio trench obviously depends on the width of the sheath.The furthermost penetration of sheath was observed in the case P = 1000 W,   = −20 V (not displayed in Figure 3).First interesting result is that an increase of substrate bias can suppress the relative deformation of sheath near the trench.The sheath deformation affects directionality of ions passing through the sheath.Therefore, it is useful to estimate the deformation of sheath before the deposition.
Ions inside the trench were our next subject of interest.The most interesting results included mean energies of ions and the difference of Ar + and Cu + behaviour.Even the case of planar substrate showed difference in energy of impinging ions.Under our range of conditions, Ar + ions lost in average from 10% to 40% of mean energy gained during their way through sheath and presheath.On the other hand, Cu + ions lost maximally 11% of gained energy.A closer examination of the situation inside the trench can explain more.As was  shown before, even in the worst case, most of the remaining electric field was located in the first quarter of the trench and still it was much lower than the field in the sheath in front of the substrate.Particles which enter the trench suddenly stop their acceleration and they should move with constant speed through the trench.Figure 4 shows the mean energies of Cu + (a) and Ar + (b) inside the trench.The conditions are   = −100 V and P = 400 W. Both species are decelerating but the effect is stronger for Ar + ions.The explanation is in the included collision processes.While only elastic collision between Cu + and Ar 0 is included, Ar + can also transfer its charge to a slow neutral Ar 0 which in fact causes significant deceleration of the ion.Ar + ions do not only lose more energy during collision but, due to higher cross section, also collide more often.The collision can change the direction of their flight and the walls catch these ions more often.This results in lower fluxes of Ar + ions to the bottom of the trench Φ Ar+ with respect to their flux to planar substrate Φ 0 Ar+ .Figure 5 shows, for example, that for   = −200 V and P = 100 W the flux of Ar + to the bottom of the trench is only approx.3% of the flux of Ar + to planar substrate.On the contrary, 35% of Cu + flux lasts down to the bottom under the same conditions.The relative flux to the bottom decreases with increasing RF power and with increasing negative bias.This corresponds to the situation of thinner sheath which is more deformed by the trench and the directionality of ions is disrupted.
Finally, we observed the fluxes of ions along the boundary of the trench.Figure 6 shows the fluxes of (a) Cu + and (b) Ar + for different biases under the condition P = 1000 W. The coordinate s follows the boundary of the trench.That is, interval (0, 2.4) represents the flat surface "to the left" from the trench; interval (2.4, 4.4) represents "left wall" of the trench; interval (4.4,4.6) is the bottom; (4.6, 6.6) is the "right wall" and (6.6, 9) is the flat surface "to the right" from the trench.Apart from the fluxes to the bottom which were discussed previously, the local maxima near the opening of the trench deserve attention.Due to high electric field, higher fluxes of ions were observed in the vicinity of the opening corners.It is interesting that the maxima are not located directly on the corner but right next to it.This effect was already observed in [17] and explained using the law of conservation of angular momentum.In addition, our observation showed that the effect of flux peaks near the trench corner is stronger in the case of thinner and hence more deformed sheaths. is a constant that describes flux of particles  to a planar substrate without trench.Coordinate  starts on the edge of the domain (0 mm), reaches the beginning of the trench (2.4 mm), then the bottom (4.4 mm), goes along the bottom (4.4-4.6 mm), then along the side (4.6-6.6 mm), and continues on remaining planar part of substrate (6.6-9mm).

Conclusions
We presented a PIC-MCC model which simulated Ar/Cu plasma behaviour in the vicinity of a two-dimensional highaspect-ratio trench.Our model provided us with detailed information about ion fluxes, number densities, energy and angular distributions, and distribution of electrostatic potential in the sheath and presheath region.The test case was IPVD of copper in argon atmosphere.A volumeaveraged model obtained the input parameters and showed the difference between Al and Cu in the magnetron chamber which was caused by different ionization energy of both species.Afterwards, the particle model showed various effects of changing RF power and bias on the sheath thickness and shape and on the amount of ions capable of reaching the bottom of the trench.Thicker sheaths, which can be achieved by higher bias or lower RF power, proved to be more suitable for depositions on the bottom of the trench.
The differences between the copper ions and ions of carrier gas were described and attributed to high cross section of charge transfer interaction between Ar + and Ar 0 .The results of our model predict the energy and angular distributions of impinging metal ions in low-pressure conditions which are characteristic for ionized plasma vapour deposition method and therefore they can be used for the optimization of this technology.

Figure 1 :
Figure 1: Scheme of the computational domain.A: source of particles, Dirichlet boundary condition  = 0 V was applied on this side.B: trenched substrate, Dirichlet boundary condition  =   .C: computational domain, period boundary condition was applied on remaining sides of the domain.

Figure 2 :
Figure 2: Validation of the global model and comparison of results for Al and Cu.(a) Electron temperature and (b) metal atoms ionization degree versus total average density of metal atoms in the chamber.

Table 1 :
A list of processes included in global model and their parameters.denotes the lowest metastable argon state and Ar * * a nearby radiative level; b A is function of electron energy.c The value is valid for neutral gas temperature 300 K and number density 2.69 × 10 19 cm −3 .d Frequency is in s −1 , electron temperature in eV, and electron number density in cm −3 . *

Table 2 :
Parameters and results of global model.Input of particle model; b free parameter; c obtained from [23]; d output of global model.