Reactive sputter deposition of TiNx films, simulated with a particle-in-cell/Monte Carlo collisions model

The physical processes in an Ar/N2 magnetron discharge used for the reactive sputter deposition of TiNx thin films were simulated with a 2d3v particle-in-cell/Monte Carlo collisions (PIC/MCC) model. Cathode currents and voltages were calculated self-consistently and compared with experiments. Also, ion fractions were calculated and validated with mass spectrometric measurements. With this PIC/MCC model, the influence of N2/Ar gas ratio on the particle densities and fluxes was investigated, taking into account the effect of the poisoned state of the target.

4 a remanent magnetic field of 13 500 G, which generates a magnetic field with a maximum radial strength of 1040 G. The dimensions of the magnets were chosen in such a way that the area of the outer magnet ring was three times higher than the area of the inner magnet cylinder. Corresponding to the classification of Window and Savvides this magnet configuration corresponds to a slightly unbalanced type II magnetron [33].
The discharge was ignited in a pure Ar plasma at 1 Pa. For a discharge current of 0.2 A, a stable discharge voltage of 260 V was measured. Without interrupting the discharge and without changing the pumping speed, an additional N 2 flow was introduced in the vacuum chamber. This N 2 gas flow was controlled with a MKS mass flow controller. Due to the stepwise increase of the N 2 flow, a corresponding change in total pressure and discharge voltage was measured. The total pressure was measured with a capacitance gauge (Baratron, Pfeiffer Vacuum), whereas the discharge voltage was logged with a normal voltmeter. Since neither the Ar flow nor the pumping speed was changed, one can assume that the change in total pressure corresponds to the N 2 partial pressure.
Ion fractions were measured with a mass spectrometer. The composition of the ion flux towards the substrate was measured with an energy-resolved mass spectrometer (EQP500, Hiden Analytical) operated in the positive-ion mode, with the orifice at ground potential. First, an energy scan at mass 36 (an isotope of Ar) was performed to determine the most probable energy of the positive ions. At this fixed energy, the amount of positive ions with a mass of 14 (N + ), 28 (N + 2 ), 36 (Ar + ) and 48 (Ti + ) was measured. Taking into account the mass abundances of these isotopes, the relative fluxes of these ions towards the orifice of the mass spectrometer could be calculated. However, due to a too strong influence of the mass spectrometer on the discharge, these measurements could not be performed at a target-orifice distance of 25 mm (as simulated). Therefore, these measurements were performed at a target-orifice distance of 70, 90, 110 and 130 mm. Since the measurements at 70 and 90 mm were comparable, one can assume that the relative fluxes at 20 mm will also be comparable with those measured at 70 mm.

Description of the magnetron in the model
The magnetron under study is based on a planar circular magnetron in a cylindrical chamber. Due to its cylindrical symmetry, the simulations can be carried out in 2d (r , z), as presented in [32]. An external resistance (R ext ) and voltage source (V ext ), which together form the external electrical circuit [32], are connected to the cathode, to create a direct current (dc). The other walls are grounded. In this manner, an electric field is generated in the discharge. In the model, the measured axisymmetric magnetic field was given as input (see section 2 above). In our simulations, we assume that the magnetron operates at room temperature (300 K), similar to the experiments. In order to study the effect of the N 2 /Ar gas ratio, the Ar partial pressure was kept constant at 1 Pa for all the calculations, whereas the partial pressure of N 2 was varied, according to the measured values (see section 2 above, and section 5.1 below).

Description of the model
As mentioned before, because of the cylindrical symmetry, the magnetron can be described in a two-dimensional space, i.e. (r , z) coordinates. However, all three velocity components must be taken into account in order to describe properly the electron gyration around the magnetic field lines, and to satisfy the energy conservation. The model developed in the present work 5 is therefore a so-called 2d3v PIC/MCC model, i.e. two dimensional in space (2d) and three dimensional in velocity (3v). The outlines of this method are given in [20], [22]- [24], [32,34], the flow chart of the model is presented in [32]. The particle movement is simulated with the PIC method (section 4.1), the collisions are treated with the MCC module (section 4.2), and plasma-surface interactions, such as sputtering, electron emission, species reflection and the effect of poisoning, are accounted for (section 4.3).

PIC method
The real particles in the discharge are represented by a limited ensemble of so-called SPs. Each SP has a weight factor W that specifies the number of real particles it represents. This is done to reduce the computation time. The program starts with an initial number of SPs and fixed number densities of the electrons, the Ar + ions, the N + ions and the N + 2 ions, i.e. 50 000 SPs for electrons, 40 000 SPs for ions and a number density of 10 15 m −3 . Initially, the SPs are uniformly distributed and are assumed to have Maxwellian velocities.
Instead of calculating the field quantities on every SP itself, they are calculated on a grid to reduce the calculation time. The SPs are first-order weighted [22] to obtain a charged particle density landscape on the (r , z) grid, with r = 0.5 mm and z = 0.1 mm.
The electrical potential V is calculated from the charge density using the Poisson equation. However, in order to simplify the solution of this equation, we use the superposition principle, so that the electrical potential can be presented as the sum of the potential only due to the space charge (V P ) and the potential only due to the cathode voltage (U 0 ): where V L is the dimensionless potential caused by an applied voltage with magnitude 1 V [35]. V P can be found as a solution of the Poisson equation, which in (r, z) coordinates reads 1 r with q being the elementary charge, 0 the dielectric constant and n i and n e the ion and electron densities, respectively. Equation (2) can be solved with zero-potential boundary conditions, i.e. V P = 0 at the surfaces z = 0 and z max and r = r max . Due to the cylindrical symmetry, the physical boundary condition at r = 0 is In [32], it is shown how the voltage drop between the electrodes (U 0 ) is calculated, resulting from the coupling of the external circuit to the plasma. In the simulation, an external circuit consisting of a constant external voltage source and resistor is applied to limit the plasma current to the desired current-voltage (I-V) regime. This external circuit is explained in [32]. V L is the solution of the Laplace equation with the applied potential boundary condition, i.e. V L = 1 at the cathode surface and V L = 0 at the grounded walls. Moreover, in the gap between the cathode and the grounded wall at z = 0, V L is assumed to decay linearly from 1 to 0 with the distance from the cathode. The Laplace equation needs to be solved only at the beginning of the simulation, which simplifies 6 and accelerates the calculation of the Poisson equation. Equations (2) and (4) are discretized using a standard five-point stencil [36] and solved on the grid. The electric field E is obtained as E = −∇V.
(5) This equation is discretized by the central finite difference method. On the boundaries, forward or backward differencing is applied.
This new electric field, together with the applied electric and magnetic fields, moves the SPs, according to Newton's equations of motion: where m is the mass, v the velocity, q the charge, x the position, E the electric field and B the magnetic field. These equations are discretized using the central finite difference method, resulting in the so-called leap-frog method [22]. The v × B rotation term is treated according to the algorithm suggested by Boris [37]. After these new positions and velocities are obtained, the procedure is repeated, until the plasma is at steady state.

MCC module
The particle collisions are treated in the MCC module. At the middle of each time step, the probability of the kth collision type is with n tar being the density of the target species and σ k (E i ) the collision cross section as a function of the energy of the incident particle E i . This probability is compared to a random number, uniformly distributed in the interval [0, 1]. If this random number is smaller, then the collision takes place, and the particles receive new positions and velocities according to the collision type. Details about their treatment for a magnetron discharge in argon are given in [23,34]. The list of the considered collisions in the Ar/N 2 discharge, and their rate constants and references to their cross sections is given in table 1. Besides elastic collisions with Ar atoms, and electron impact ionization and excitation of Ar (ground state and metastable atoms), also electron impact ionization, excitation (to four different excited levels), dissociative ionization and dissociation of the N 2 gas molecules are included. Ionization of N is not considered, because N has a lower density, and hence its ionization is less important. Also, the density of N + is low, so recombination of N + is also not included. Elastic collisions of electrons with N 2 are also omitted, due to the lower density of N 2 , so the momentum change of the electrons by elastic collisions with N 2 is small compared with Ar. As far as the heavy particle collisions are concerned, elastic scattering of Ti + with N 2 and N is included, because, in constrast to electrons, a considerable amount of energy is transferred. Charge transfer of Ar + with N 2 is included, due to its very high rate constant. The model also contains elastic and charge transfer N( 2 ) + collisions with Ar, N 2 and N.
Note that not only electrons, Ar + ions, N + ions and N + 2 ions, which are assumed to be present from the beginning, are included in the model. Also fast Ar atoms (that originate from elastic collisions, including symmetric charge transfer collisions, with Ar + ions [38], i.e. reactions (16) and (17) from table 1), metastable Ar atoms, fast Ti atoms, Ti + ions and fast N atoms, created out of the initial species (by plasma reactions or by sputtering) are considered. Slow Ti and N atoms are described in the model with balance equations (see section 4.3). Some cross sections of ion-neutral collisions are described with the Langevin-Hasse model [39,40], assigned 'LH' in table 1: where α p is the polarizability, e the electron charge, 0 the dielectric constant of vacuum, µ the reduced mass, g = |v i − v n | the relative precollision velocity, with v i and v n the ion and neutral velocities, respectively, and β ∞ is the value of the dimensionless impact parameter β, for which the deflection angle is negligibly small [41]. This value is set to 3 for Ar, N 2 and N [40]. The polarizability for Ar is 11.08a 3 0 [42], 18.24a 3 0 for N 2 [43] and 7.5a 3 0 for N [42], where a 0 is the Bohr radius.

Plasma-surface interactions, implemented in the model
Since magnetron discharges are mainly applied as sputter devices, the sputtering of the target is included in the model. This is done by calculating the sputter yield Y , each time an ion or atom hits the cathode surface. The sputter yield for an incident particle i with energy E i is described by the empirical formula of Matsunami [44] Y where U s is the sublimation energy of the cathode, E th the threshold energy and the other parameters are properties of the cathode material, as described in [44].

9
The sputtered Ti atoms (and N, in the case of a poisoned target, see section 4.4) are followed with the PIC/MCC method until they are thermalized in order to reduce the computation time. However, the Ti (and N) atom density is important in picturing the deposition process. Therefore a compromise is found between computational effort and accuracy in treating the thermalized Ti (and N) atoms as a fluid, and hence their density n slow (r, z) is calculated with the diffusion equation where D is the diffusion coefficient of Ti or N atoms in Ar gas, and r prod and r loss are the production and loss rates, respectively, of the Ti or N atoms. The production mechanisms of slow Ti atoms are defined by transfer from the fast (sputtered) Ti atoms, i.e. thermalization by elastic collisions with Ar (reaction (31) from table 1) and thermalization at the walls, whereas the loss mechanisms of slow Ti atoms are electron impact ionization (reaction (7) from table 1), and Penning ionization (reaction (23) from table 1). Slow N atoms are produced by thermalization of fast N atoms at the walls, as well as by thermalization of these fast N atoms by elastic collisions with Ti + ions (reaction (30)), N + ions (reaction (37)) and N + 2 ions (reaction (43)). Other production mechanisms are dissociative ionization (reaction (13)) and dissociation of N 2 gas molecules (reaction (14)), and dissociative recombination of N + 2 (reaction (15)). Chargetransfer reactions of N + ions can also cause the creation of slow N atoms, i.e. with Ar gas atoms (reaction (34)), N 2 gas molecules (reaction (36)), N atoms (reaction (38)). Slow N atoms can be lost by elastic collisions with N + ions, including charge transfer, and by charge transfer with N + 2 ions (reaction (44)), because these processes give rise to the creation of fast N atoms.
The diffusion coefficient (in m 2 s −1 ) is calculated from the rigid sphere approximation [45] where T is the gas temperature (K), M 1 and M 2 are the masses of Ti (or N), and Ar (g mole −1 ), p is the pressure (atm), and d 12 is the collision diameter (10 −10 m), given by (d 1 + d 2 )/2. The collision diameter of Ti is 2.684 × 10 −10 m [46], 3.298 × 10 −10 m for N [47], and 3.542 × 10 −10 m for Ar [47]. The second plasma-surface interaction taken into account is ion-or atom-induced secondary electron emission, characterized by the secondary electron emission coefficient (SEEC). The SEEC describes the number of secondary electrons produced by an atom or ion hitting the cathode surface. Note that secondary electron emission at the other walls is less important, so it is not included. Unfortunately, a wide range of different SEEC values are reported for ion bombardement of Ti (from 0.075 [48] to 0.148 [49]). Moreover, some of these values describe the effective secondary electron yield, which is dependent on both SEEC and reflection coefficient (RC) [50] (see below). Therefore, the exact SEEC is not known.
Another plasma-surface interaction regards an electron striking a wall, after which it can be reflected or adsorbed, characterized by the RC. In magnetrons, the electron density near the cathode is much higher than at the other walls, so only interaction of electrons with the cathode surface is important. The RC is, however, hard to measure and to our knowledge no values are reported in the literature. Since both the RC and the SEEC directly influence the cathode current and voltage [51], they are slightly adapted in the model, so that the calculated currents and voltages can be compared with experimental values. Note, however, that these coefficients were kept the same for the different I-V combinations investigated for sputtering of Ti in pure Ar.
Similarly to electrons, heavy particles can also be reflected or adsorbed when hitting a wall, depending on the sticking coefficient (SC). The SC of a reactive N atom is mostly assumed to be 1 [26], whereas the SC of Ti is found to be dependent on target-substrate distance [52]. In accordance to these reported values [52], an SC Ti of 0.5 is chosen in our model. The SC of the background gases is assumed in the model to be 0. For Ar being an inert gas, a zero sticking approximation is justified. The zero sticking for N 2 , on the other hand, is chosen based on the following: the partial pressure of N 2 is used in the code as measured by experiment. In the simulation, this constant partial pressure corresponds to the assumption of a zero effective SC for N 2 . Moreover, the SC of N 2 is very low, and its influence on the calculated N 2 density is therefore negligible, which also justifies the assumption.
When a heavy particle is reflected, a distinction is made between elastic and inelastic reflections depending on the thermal accommodation coefficient, defined as where E i is the mean energy of the incoming particles, E r the mean energy of the reflected particles and E w the mean energy of the reflected particles in thermal equilibrium with the wall (2kT w ). This means that if E r = E w , all the energy of the incoming particle is transfered to the wall, and an inelastic reflection occurred. In this case, α is equal to 1. If, however, E r = E i then the energy is conserved, and an elastic reflection occurred. In this case α is 0. In the simulation, α is assumed to be 0.5 [53], and a random number is generated. When the random number is lower than α, an inelastic reflection occurs, and if the random number is higher than α, an elastic reflection occurs.

Effect of poisoning on the plasma-surface interactions and on their implementation in the model
When a reactive gas, like nitrogen, is added to the argon background gas, N( 2 ) + ions can be implanted in the Ti target, nitrogen molecules and atoms can be chemisorbed at the target surface, and chemisorbed species can be knock-on implanted into the target. Subsequently, they react with the target atoms, to form a TiN x layer, similar to the mechanism in oxygen [30]. This surface modification process is called 'poisoning', and influences the plasma-target interactions (i.e. sputtering and secondary electron emission), and therefore all of the plasma properties. The transition from the so-called metallic to poisoned condition happens via a hysteresis [30,31]. In the present work, we used N 2 gas flows for which the target is completely in poisoned mode, in order to avoid (i) the simulation of the hysteresis and (ii) having to deal with a partially reacted Ti target. Both Ti and N atoms can be sputtered from a poisoned TiN x target, but their sputter yield is lower than for a metallic Ti target. According to the values reported in [26], the sputter yield of Ti from a fully poisoned target is lowered with a factor of 6.4 compared with the sputter yield of Ti from a metallic target. On the other hand, the sputter yield of N from a fully poisoned target is 4 times higher, compared to the sputter yield of Ti from a fully poisoned target. These changes in sputter yields for the poisoned target apply to all of the bombarding species (i.e. Ar + , Ar f , N + 2 , N + and Ti + ; see later, section 5.4 and tables 2 and 3). Moreover, due to the target surface modification the SEEC is altered that results in a target voltage change [50]. In the model, the poisoning of the target is therefore described by changing the SEEC value, depending on the N 2 gas flow. The SEEC value is adjusted, in accordance with the range reported in [50] (see values below).

Input values and calculated I-V characteristics
In our calculations, an RC of 0.1 and an SEEC of 0.075 yielded calculated I-V values in good agreement with the experimental data in pure Ar. Starting from these coefficients, an Ar/N 2 mixture is evaluated in a constant current regime of 0.2 A (accomplished with V ext = −600 V and R ext = 1500 ). The applied magnetic field has a maximum radial magnetic field strength of 1040 G. To study the effect of the N 2 /Ar gas proportion, the Ar partial pressure is kept constant at 1 Pa for all calculations, whereas the N 2 partial pressure is increased as 0.03, 0.06, 0.13, 0.19 and 0.26 Pa. Under these conditions, the target is always fully poisoned, so that the hysteresis behaviour can be avoided. To summarize, only the partial pressures of the Ar and N 2 gases, V ext and R ext , the magnetic field, the SC, RC and SEEC (discussed in sections 4.3 and 4.4), and the cross sections and rate coefficients of the various collisions in the plasma (see table 1) are input values in the model. All the other plasma characteristics, including the results presented below, are calculated self-consistently.
The SEEC alters as a consequence of poisoning, and in the case of a TiN x target, the SEEC decreases [50]. Note that the SEEC values applied in the model comprise the SEEC values of all different incident species, to avoid complicating the model with different uncertain parameters. From figure 1, it is clear that the overall SEEC is lowered in the model with increasing the N 2 pressure. However, in reality, the SEEC values of individual species will probably decrease first, but then remain constant once the poisoning is complete. Nevertheless, the proportion N( 2 ) + /Ar + will increase with N 2 flow, and because the SEEC of N( 2 ) + is much lower [54] than the SEEC of Ar + the overall SEEC will indeed decrease with increasing N 2 flow.
The SEEC values used, and the calculated currents and voltages as a function of N 2 partial pressure are illustrated in figure 1, along with the experimental values.

Calculated electron, ion and atom densities
The externally applied magnetic field traps the electrons in an area close to the cathode. Most electrons are trapped in the region where the radial magnetic field is at a maximum, causing a peak in the electron density. This is clear from figure 2, presenting the electron density in the simulation area (r, z), for the case of 0.26 Pa N 2 . Similar profiles for the other N 2 partial pressures were obtained, but they are not shown. Indeed, since the SEEC only varies within a limited range due to gas composition (see section 5.1), it will not have a large effect on the electron density.
The localized electrons ionize neutrals, leading to similar density profiles for the Ar + ions, N + ions, N + 2 ions and Ti + ions, as shown in figure 3, for an N 2 partial pressure of 0.26 Pa. It is clear from this figure that the Ar + ions are the dominant positive ions, with a maximum density of 7 × 10 17 m −3 , which is only slightly lower than the maximum electron density (see figure 2). The N + 2 ions reach a maximum density of about 1.2 × 10 17 m −3 , which is a factor of almost 6 lower than the Ar + density, despite the fact that the N 2 partial pressure is only a factor of 4 lower than the Ar partial pressure. This is attributed to the fact that the N + 2 ions are lost more efficiently (by dissociative recombination with electrons) than the Ar + ions. The N + ion density is still more than an order of magnitude lower (with a maximum density of about 9 × 10 15 m −3 ), which can be explained by the rather low dissociation degree of N 2 (i.e. the N atom density is also much lower than the N 2 density, as will be shown below). The Ti + ions have an even lower density, with a maximum of only 1.8 × 10 15 m −3 , because these species do not originate from the background gases, but only from ionization of the sputtered atoms. As most of the fast Ar atoms originate from charge transfer reactions of Ar + ions, the typical Ar + ion peak profile appears for the fast Ar atoms as well, see figure 4. Its density is quite high, compared to the ion densities, but it is still two orders of magnitude lower than the overall background Ar gas density, which is of about 2.4 × 10 20 m −3 . Ti atoms originate from sputtering the cathode target, and therefore, the Ti density has a maximum near the cathode. Its overall density is 3-4 orders of magnitude lower than the Ar atom density. N atoms are also sputtered from the poisoned target, but most N atoms are created in the plasma instead of by sputtering (by reactions (13)-(15), (34) and (36) from table 1). As a consequence, the peak near the cathode is not so pronounced as in the case of the sputtered Ti atoms, and the N density is characterized by a broad profile throughout the discharge. Also, the N density decreases towards the walls, because a SC equal 1 was assumed for N. Because N is not only created by sputtering but also by plasma reactions, the maximum value is almost an order of magnitude higher than the sputtered Ti atom density. Comparing the average N atom density, which is about 1.4 × 10 17 m −3 , to the N 2 molecule density, which is about 6.2 × 10 19 m −3 , reveals that the dissociation degree is in the order of 0.2%.
To investigate the influence of the N 2 /Ar gas proportion on the various plasma species, 1d density profiles on a line perpendicular to the cathode at the peak density (i.e. r = 13.5 mm)  Figure 5 illustrates the results for the Ti and N atoms. As a consequence of the difference in the Ti sputter yield for a metallic or a poisoned target (as mentioned in section 4.4, and as will be shown below in section 5.4), the Ti density drops significantly upon N 2 addition. However, when adding more N 2 gas, the Ti density remains constant. The N density rises  logically with N 2 gas amount. The explanation for the profiles of the Ti and N densities was given in the paragraph above. The influence of the N 2 /Ar gas proportion on the 1d ion density profiles is shown in figure 6. When raising the N 2 flow, the N + and N + 2 densities increase and the Ar + density decreases. This is explained as follows: the charge transfer reaction of Ar + with N 2 , which causes the production of N + 2 ions and the loss of Ar + ions, has a high rate constant (see reaction (21) of table 1). As a consequence, an increased N 2 gas amount leads to a higher N + density, but to a lower Ar + density. The Ti + density drops approximately a factor of 10 when N 2 is added, similar to the Ti atom density (figure 5), caused by target poisoning.
Mass spectrometric measurements were carried out to determine the ion fractions at 7 cm from the cathode, for different N 2 gas concentrations. However, in order to save computational time, the size of the simulated magnetron reactor was limited to 2.4 cm. Hence, the ion fractions are calculated at 2 cm from the cathode. Nevertheless, it is assumed that the ion proportions will not vary much in the bulk of the plasma. From the measured and calculated ion fractions, presented in figure 7, we conclude that a good agreement with experiment is found, and similar trends as for the ion densities are found. Only for the Ti + ions is the agreement less satisfactory.

Calculated ion and atom fluxes to the cathode
The ions accelerate towards the cathode by the applied electric field, where they can sputter the target. Also some neutrals contribute to the sputtering (see below). In figure 8, particle fluxes to the cathode target, as a function of radial position, are presented for different N 2 concentrations. The peak at 13.5 mm is a direct consequence of the maximum in the density profiles (see section 5.2 and figures 3 and 4). When the N 2 partial pressure is increased, the N + 2 and N + fluxes increase and the Ar + flux decreases, in analogy to their densities (see section 5.2 and figure 6). The Ti + flux drops approximately a factor of 15 when N 2 is added. As was the case for the Ti + density (see section 5.2), this is a consequence of the lower Ti sputter yield when the target is poisoned (as will be shown in section 5.4). The fast Ar atom flux to the target is rather independent of the N 2 amount.

Calculated ion and atom contributions to sputtering and total sputter fluxes
Depending on the magnitude of the fluxes, the corresponding energies, masses and atom numbers, the above-mentioned particles contribute to the sputtering of the target (see equation (9) above). In tables 2 and 3, the relative amount (in %) of Ti and N sputtering, respectively, created by each of these ions and atoms, is summarized.
For both Ti and N, most sputtering is caused by bombarding Ar + ions, followed by Ar f atoms. As mentioned before, these Ar f atoms originate from elastic collisions (reaction (16)), including symmetric charge transfer collisions (reaction (17)), with Ar + ions [38]. The N + 2 ions only play a role at high N 2 partial pressures, and more for sputtering N than Ti, because of the smaller mass differences, and hence the higher sputter yield. The role of N + and Ti + ions can be neglected under the investigated conditions, with a contribution of about 0.1-0.5%. The order  in contribution is a consequence of the magnitude of the fluxes, as seen in figure 8. In general, raising the N 2 partial pressure causes an increase of the N + 2 and N + contributions and a decrease of the Ar + contribution to the sputtering of both Ti and N. This can also be explained by the dependance of the N + 2 , N + and Ar + fluxes on the N 2 partial pressure, as seen in section 5.3. From table 2, it is also clear that, when no N 2 gas is present, the contribution of Ti + to sputtering the metallic Ti target is a factor of 15 higher. This is a direct consequence of the Ti + flux (section 5.3 and figure 8).
The total fluxes of sputtered Ti and N atoms are shown in figure 9. The maxima of the ion and atom fluxes towards the cathode cause a maximum in the sputter flux. This localized erosion of the target creates the so-called 'race track'. When no N 2 gas is added, the sputtered Ti flux is approximately a factor of 8 higher than after adding N 2 , and hence poisoning the target, causing the sputter yield to decrease. Once the target is poisoned, the sputtered Ti flux remains constant and does not drop further upon N 2 addition. The sputtered N flux is higher than the sputtered Ti flux due to the higher sputter yield of N (see above, section 4.4). Moreover, in contrast to the Ti flux, the sputtered N flux increases slightly with increasing the N 2 partial pressure. This is a consequence of the dependance of the N or Ti sputter yield on the different incoming species. Indeed, sputtering N with N + 2 or N + ions has a larger yield than with Ar + ions, whereas sputtering Ti is less dependent on the bombarding ion type. With increasing the N 2 flow, the N + 2 and N + fluxes increase (see section 5.3 and figure 8), causing an enhanced sputtering of N. Calculated sputtered Ti (dashed lines) and N (solid lines) fluxes from the cathode for different N 2 partial pressures as indicated in the legend and at an Ar partial pressure of 1.0 Pa. The lower part of the Y-axis is stretched for clarity, because the sputtered Ti flux in pure Ar gas is much higher than in an Ar/N 2 mixture. Note that the Y-axis (r = 0) corresponds to the symmetry axis of the cylindrically symmetrical reactor.

Calculated Ti and N fluxes to the substrate
In order to obtain a better insight in the deposition of TiN x films, the calculated fluxes of Ti and N atoms to the substrate are presented in figure 10.
The Ti flux to the anode drops when N 2 is added, and remains constant afterwards, in correspondance to the sputtered Ti flux (see section 5.4 and figure 9). The Ti flux to the anode is characterized by a similar radial peak profile as the sputtered Ti flux, implying that the deposited Ti in the film will be non-uniform. The broadening of the peak profile is a consequence of diffusion of the sputtered Ti atoms through the plasma. In analogy to the sputtered fluxes (see section 5.4 and figure 9), the N flux to the anode is higher than the Ti flux almost everywhere. The N flux has lost its radial peak profile, as was also clear from the N density (see section 5.2 and figure 4), because more N is created by plasma reactions than by sputtering. The N flux increases with N 2 amount, similar to the N density (see figure 5).
When the absolute values of the fluxes of Ti and N to the substrate are compared (figure 10), and taking into account the assumed SCs of 0.5 and 1, respectively, most of the deposited TiN x film would have a stoichiometry x much greater than 1. This implies that the assumed SC of N atoms must be much lower, probably of the order of 0.1, as reported in [55]. On the other hand, as stated in [26], the effective SC of a species will be a function of the coverage of the different species on the deposited film. Indeed, with increasing N coverage, the effective SC N will be lower. Therefore, it would be advisable to couple our PIC/MCC model to a surface model, such as [27,56], to account for this lower effective SC of N. Figure 10. Calculated Ti (dashed lines) and N (solid lines) fluxes to the anode for different N 2 partial pressures as indicated in the legend and at an Ar partial pressure of 1.0 Pa. The lower part of the Y-axis is stretched for clarity, because the Ti flux to the anode in pure Ar gas is much higher than in an Ar/N 2 mixture. Note that the Y-axis (r = 0) corresponds to the symmetry axis of the cylindrically symmetrical reactor.

Conclusion
To simulate the physical processes in a magnetron discharge during reactive sputter deposition, a 2d3v PIC/MCC modelling approach was applied. With this model, we are able to calculate in detail the densities and fluxes of the various species in the whole reactor. The model was validated by comparing calculated ion fractions to mass spectrometric measurements. Although this modelling approach is very accurate and self-consistent, it suffers from a long calculation time. Consequently, it is not so suitable for large reactors. However, its precise results can be used to understand certain mechanisms, and hence to improve simpler, but faster models, such as fluid and analytical models, which can deal with larger reactors.
This model accounts for poisoning of the titanium target when N 2 gas is added, using N 2 pressures for which the target is completely poisoned. However, the poisoning process of the target at low N 2 pressures, i.e. the hysteresis behaviour, is not described in the model. Reaching steady state at low N 2 pressures takes several seconds to minutes, and therefore a PIC/MCC approach will not be suitable.
With respect to the deposition of TiN x films, a stoichiometry x much larger than one is predicted, when SCs of Ti and N of 0.5 and 1, respectively, are used. Therefore, the effective SC of N will probably be much lower than 1, and will be a function of the film coverage. In the future, we plan to couple the PIC/MCC model to a surface model to account for this. Also, the sticking of molecular N 2 gas will be included in the surface model, because N 2 can also contribute to the deposition of the TiN x film due to its high density. Nevertheless, this model is already able to provide accurate information on the processes that occur in the plasma. It can be used as an extension or even validation for experiments and simple models, to get a step closer to a complete view of reactive magnetron sputter deposition.
In the future, this PIC/MCC model can be extended to other reactive gases (e.g. Ar/O 2 ) or combinations (e.g. Ar/N 2 /O 2 ), and to other electric (e.g. rf, pulsed dc, maybe HIPIMS) or magnetic (e.g. unbalanced) configurations.