Multiscale Molecular Dynamics Simulation of Plasma Processing: Application to Plasma Sputtering

Molecular dynamics is an atomistic tool that is able to treat dynamics of atom/molecules/cluster assemblies mainly in the condensed and liquid phases. The goal of the present article is to provide a new methodology for describing all phenomena of plasma processing and beyond such as gas phase chemistry as well. Simulations of condensed matter and liquid processes by molecular dynamics are now readily accessible provided the interaction potentials are available, so quantitative parameters can be deduced as diffusion coefﬁcient, etc. The situation is less clear for gas phase processes while they operate on larger space and time scales than for condensed phases and at lower specie densities. The present article is proposing a new methodology for describing plasma core interactions in taking into account experimental space and time scales. This is illustrated on a plasma sputtering process and deposition in a single simulation.


INTRODUCTION
Multiscale modeling of plasma processing and more generally of chemical processing can be considered as a "holy Grail" quest.Numerous work have been carried out since the 90' .The ultimate goal is to treat in a "quick way" all the length and time scales of the processing, i.e., from nanometer to centimeter or even meter and from picosecond to hours.It is worth noting that such approaches have mainly concerned interaction with materials.Coupling between scales is achieved in a hierarchical way: The fundamental growth mechanisms, both in gas phase and at the materials interface, occur on the atomic scale.But, the geometry of the deposition/etching/treatment reactor and the associated operating parameters directly affect the chemical composition of the gas/plasma and the temperature at the growth surface, if any.The properties are, in turn, controlled by both atomic-and microstructural-scale features.This approach provides the missing link between chemical vapor deposition reactor design/operating conditions and the material structure/properties in the case of film growth.A tricky challenge is to use atomic scale simulation methods for enabling multiscale analysis.Kinetic Monte-Carlo methods are now relevant for describing a wide range of time and space scales [9,12,16,[22][23][24][25], including plasma materials processing reactors [26][27][28].Recently Molecular Dynamics simulations (MDs), which are calculating all the trajectories of an ensemble of particle, have overcome the barrier of describing long time phenomena or rare infrequent events [29][30][31][32][33][34].Indeed, MDs are highly demanding computational resources while "exactly" solving Newton or Langevin equations of motion, which, a priori prevents to reach long time dynamics without including additional ad-hoc powerful procedures.Moreover, the recent availability of (reactive) many-body potential [35,36] has increased the interest of carrying MD simulations of complex systems both in materials science and in plasma/gas phases.
For going further, the present manuscript aims to provide a method for a full description of all steps of plasma processing, from gas phase to thin film growth or materials treatments.

Molecular Dynamics Principles
Molecular Dynamics simulation (MD) is a powerful techniques for describing atomic and molecular processes at the nanoscales.It is solving the equation of motion using Newton equations of motion (or Langevin dynamics which includes energy dissipation through an additional friction term): where − → r i (t) is the position of atom i at time t with mass m i , and V is the interaction potential between all N involved species.So solving such a set of equations requires the knowledge of the interaction potentials V and a set of initial conditions, preferably matching experimental situations.It is widely recognized, using massively parallel programming and computers, that MDs are able to solve multi-million atom system dynamics (which corresponds to a volume of the order of 100 × 100 × 100 nm 3 ) over several nanosecond time-scale with many body-interactions.In such a form MDs, today, are not able to go further in size and time scale.Recent advances either by coupling Monte-Carlo approaches or accelerated dynamics allowed to move time scales up to several seconds [34,37].Nevertheless, for addressing large reactor scales, it is necessary to find a way for rescaling experimental conditions toward simulation sizes and accessible computational time.
First of all, it is interesting to point out in which situations MDs are able to give relevant macroscopic information, ideally only depending on the interaction potential quality.
In condensed and liquid matter, distances between atoms or molecules are small enough for defining a simulation box identical to a representative part of the real system.For illustrating this, let us consider a single crystal, e.g., platinum.Pt density is 65 at.nm −3 so in a 10 × 10 × 10 nm 3 cell there are already enough atoms for calculating average macroscopic quantity with enough precision, provided accurate interaction potentials are available, which is the case for Pt.For example self Pt diffusion coefficient can be calculated as well as numerous other quantities.In liquid, say water, there are 33 H 2 O molecules per nm −3 , so it is also possible to have a large number of molecules in a simulation box for reaching macroscopic quantities with enough precision, also if interaction potentials are of sufficient precision.In this case, time -dependent quantities are also deduced with enough accuracy while no hypothesis on space and time scales are done.
The situation in diluted media, among them plasmas, is quite different.For example, the density of a gas at atmospheric pressure, P = 1. 10 5 Pa, is 2.4 10 −2 nm −3 , so it is quite difficult to have a tractable simulation box with a significant atom/molecule number, say 50-1,00,000.Moreover low pressure plasma processing involves pressure as low as 0.1-10 Pa, which drastically lowers the possibility of handling molecular dynamics.
A second difficulty is arising when considering plasmamaterials interaction: the fluxes of species on material surfaces.Typical fluxes for deposition/etching/surface treatment are typically of the order of 1. 10 15 cm −2 .s−1 , which is close to one surface layer of a flat material surface per second.Such low flux is not reachable by MDs: it is not possible to directly wait for 1 s when the integration time interval for solving Equation ( 1) is in the range 0.1-5 femtosecond.This problem can be simply overcome by considering that the time interval for releasing an atom/molecule toward a surface should be large enough for allowing energy relaxation between two consecutive impacts on the surface.When choosing correctly this time interval, the results should not change when increasing it [38,39].But the true time of experiments is not recovered, and thus consequently macroscopic rates as etching, depositions rates, etc., also do not.

Pushing MDS Toward Multiscale Description of Plasma Processing
Considering all points listed above, realizing multiscale modeling based on MDs requires to be able to describe gas phase reactions in a realistic way and with a suitable time scale.
First of all, plasma is a complex medium including electrons, ions, ionic or neutral radicals, and excited states (electronic, rotational, and vibrational) of these latter species.The most difficult species to include in MDs are electrons and excited states [31].Nevertheless, interaction potentials aiming at describing electron motion in MDs have been proposed: these are eFF [40,41] and eReaxFF [42] force fields.Both are suitable for including electrons in simulations.eFF allows excited atomic states to be treated so monoatomic plasmas are reachable, while eReaxFF is explicitly treating electrons in chemical reactions using the reactive potential.At this stage electronic excited states can be included in MDs, and chemical reactions involving (ionic) radical, electron collisions, etc, are expected to work.A pending question is the description of vibrationally excited molecules and radicals.Vibrational excited states are quantified and their chemical reactivity, as occurring in plasma chemistry, is thus difficult to capture using classical MDs.Nevertheless some attempts have been done with a relative success [43] and are opening the way for new studies.Two ways are explored: increasing potential energy of the molecule by stretching the bonds to an appropriate length or increasing velocities of the atoms of the molecule which results in increasing kinetic energy of the atoms of the molecules resulting in oscillations with larger bond length.As a summary, Table 1 is displaying capabilities of MDs in including the nature and state of plasma species.
If we consider that all these issues specific to plasmas are solved or are giving a first approximation for starting MDs of physics/chemistry inside a reactor, a last work should be done: how to describe chemical reactions and or collisions at low pressure for which distances between interacting species are large, even when mean free path is lower than reactor size for allowing collisions and then reactions.If we recall that during two collisions, neutral species are traveling in a straight line while ions have an ambipolar diffusion in the electric fields (such a motion can be simulated using MDs), nothing happens except radiative de-excitation, if any.So the way for having a correct description of the plasma core is to build a simulation box with the same collision number and with a distance between plasma species greater than twice the largest interaction cutoff length.Recall that, the collision number is proportional in experiments and simulations, to the pressure-distance product through the relation: T g is gas temperature, n is the number of collisions that take place in the gas, d is the reactor or simulation box size, P is the gas pressure in the is the reactor or simulation box, and σ See text for explanation of the different factors.Reprinted from Voter [31].
is the collision cross section.Equating the collision number in simulation box (SB, label sim) and in the reactor (label exp) leads to the relation: Noticing that the SB volume can written as V sim = S sim d sim where S sim is the area normal to d sim .If d sim is along z axis then S sim is in the xy plane.Using P sim = N sim V sim • k B T g , with N sim being the number of atoms in the SB, Equation (3) can be rewritten as: which does not depend on d sim .Nevertheless, d sim can be calculated keeping in mind that the species in the gas/plasma phase should not be closer than the largest interaction cutoff r cut used in the simulations.The distance l between gas/plasma phase species can thus be estimated to: > r cut , and thus And so the correct number of species in the simulation box can be deduced for matching experimental conditions.Moreover, if concentration, temperature, etc. gradients are present in the reactor, the gas/plasma can be divided in representative cells as in Figure 1, including substrate for thin film deposition/etching/treatments, if any.Each of these cells is associated with a simulation box of a size preserving the collision number, with different parameters, pressure and/or temperature for example.The dynamics is calculated in the considered box and at a relevant simulation time, the created species are allowed to move toward a neighboring simulation box describing the neighboring cell in the experimental setup (see Figure 1).So the full simulation box will be made of individual boxes communicating in the same way as in the experimental reactor.Each individual boxes is the "model" of its experimental counterpart and is appearing in the same ordering as the cells in the experiments.This mapping is displayed as a simplified 2D picture in Figure 1.
When succeeding in this mapping, there is a remaining issue to be solved: converting the simulation time to the experiment time, in the gas/plasma phase.
It should be recalled that during the flight between two collisions, nothing happens that changes velocities and direction of the particles and that the velocity in simulations and experiments should be the same.As shown in Figure 2, if we are always considering that the collision number is the same in an experimental cell and its corresponding simulation cell, the unaffected velocity v between two collisions can be written as: where we define the traveled distances d exp , d sim and elapsed times t exp , t sim between two collisions, respectively for experiments and simulations.So d exp , d sim are the collision (interaction) mean free paths.ρ is the interaction distance, assumed for simplicity to be the same for the plasma species.Generalizing to different interaction lengths ρ i is straight forward.Thus the expected experimental time can be deduced from MDs as: When particle are interacting, then the interaction duration are the same in both experiments and simulations, because distances between species become identical: Accessing the "experimental time" will thus allow to directly calculate the reaction rate coefficients in the gas/plasma phase, while usually classical MDs is rather able to calculate reaction probability, which can be connected to reaction rates using relevant approximations.This description of the plasma core reactivity is thus providing, in a space and time resolved way, the consistent source of species that will impact on the substrate surface and reactor walls leading to deposition/etching and/or other relevant surface interactions.
Such methodology is not limited to plasma processing, and it can be applied for any chemical processing, provided interaction potentials are available.

APPLICATION TO PLASMA SPUTTERING
For validating the method on a simple case, this MD procedure is applied for describing plasma sputtering deposition, including sputtering, transport of the sputtered atoms through the plasma phase and deposition on a substrate.At this stage including all plasma processes (electron collision using electric force fields, excited states, etc.) is beyond the scope of this article.The sputtering of a platinum target in an argon plasma with subsequent deposition onto a porous carbon is addressed, while both processes have been previously and separately studied using MD [44,45] and experimentally [46,47].
The simulation box is thus built as a porous carbon substrate, a platinum target, an argon gas in between.The simulations box size is 12 × 12 × 136 nm 3 , with periodic boundary conditions in the xy plane.The total number of atoms in the SB is 170,056, including 68,160 carbon atoms, 96,100 Pt atoms and 5,796 argon atoms among them 500 are randomly selected during the simulation to become Ar ions directed toward the Pt substrate with the cathode fall kinetic energy.The total Ar atom number is calculated using Equation (4) with experimental argon pressure P Ar = 2.5 Pa (n coll = 8), with a target to substrate distance of d = 6.5 cm [46,47] and temperature T of 300 K.
The impinging Ar ion energy on the Pt target is chosen to be 500 eV: this allows significant sputtering of Pt atoms [45].They are released toward porous carbon substrate every 2 ps.Ar gas atoms are thermostated at 300 K: this is mimicking energy transfers to the reactor wall after Ar + -Ar and Pt -Ar collisions.Pt target is divided in 3 zones.The top of the crystal is made of 0.6 nm thickness zone with motionless atoms for preventing Pt target motion under Ar ion impact.
A 7 nm thickness zone is thermostated for dissipating the energy released by the impinging Ar ions.The remaining zone, 2.4 nm thickness, is composed of free moving atoms.They are not thermostated for allowing correct escape kinetic energy of sputtered Pt [45].The model porous carbon [48] is divided in two zones: the bottom zone with motionless carbon atoms, 0.6 nm thickness prevents porous carbon substrate motion under impact of depositing Pt atoms.The remaining part is thermostated at 300 K mimicking substrate temperature control, here at 300 K.As a summary, Figure 4 illustrates the MD simulation initial configuration.
The Pt interaction potential is the Embedded Atom Model (EAM) function [45,49,50], the Tersoff potential is used for C-C interactions [51].Repulsive interaction of Ar with Pt, C, Ar are of the Moliere form Graves and Brault [52].Pt-C interaction are modeled by an improved Lennard-Jones potential [44,53].The equations of motion (1) are integrated using velocity-Verlet algorithm, with a timestep of 0.1 fs.Such a short timestep is necessary for describing Ar ion-Pt interactions.The LAMMPS software [54] 1 is used for running the simulations.Damping time for Berendsen thermostats used for dissipating energy in the plasma region through Ar gas and in the carbon substrate and platinum target is 1 ps (see Figure 3).The total calculation time is 2 ns: during the first ns, the Ar ions are impinging the Pt target and sputtered Pt atoms are traveling toward the surface.An additional calculation time of 1 ns is needed for allowing last created Pt atoms to interact with the carbon substrate.
Figure 4 is presenting the snapshots at 2 ns elapsed time (2 10 7 iterations) for this present simulation of the full sputtering deposition process (Figures 4a-c) and with the snapshots obtained using simulations described in Xie et al. [44] ( Figures 4d-f).In the latter case, the deposition occurs without Ar gas and for which the initial velocity distribution is calculated from the Thompson distribution at the substrate position [44], which is expected to take into account the collision between traveling Pt atoms and Ar background between target and substrate.The comparison between these two sets of snapshots (a-c and d-f) shows that cluster size distribution and density are the same, as well as the depth profiles.So it can be concluded that the reactor model is able to reproduce the deposition process.Moreover, it should be recalled that the Pt clusters morphologies were in agreement with experimental findings [44].
So it can be deduced that the sputtering model of sputtering reactor presented here is accounting for a correct description of sputtering, transport and deposition of Pt onto porous carbon.It should be mentioned this is showing that the transport from the target to the substrate is correctly taken into account.This is establishing that the assumption of considering the collision number as a prescribed input of the model is correct.Even if the present model does not includes electron collisions, this is a first step that is assessing the interest of the approach.When there is no electron, say in chemical vapor deposition, this work gives the way for treating reactivity in the plasma core and the interaction with the materials.This will be achieved by a simulation box, with relevant dimensions calculated using Equations (4-5) filled with molecules of interest, interacting with the chosen substrate.

CONCLUSION
MDs capabilities in materials science are increasing with the availability of new reactive many-body potentials, macroscopic parameters are readily accessible and comparison with experiments is more and more efficient.For gases and plasmas, the connection between MDs and experiments is difficult due to the too diluted matter preventing direct realistic large scale simulation to be carried out.In the present article, when highlighting that considering the collision number in both MDS and experiments, a direct link can be established between MDS and experimental results allowing recovery of macroscopic time dependent coefficients (reactions, deposition, etching rates, etc.).A general procedure for describing plasma (chemical) processing is presented, suggesting a possible way for recovering reaction rates based on distance scaling between experiment and simulation conditions.As example, a sputtering model is presented which includes sputtering, transport, and deposition making use of the conservation of the collision number.It correctly describes the overall sputtering process, consistent with previous experimental findings.Future work will include electron collisions in the gas phase.

FIGURE 1 |
FIGURE 1 | 2D schematics of the mapping of (A) a plasma reactor illustrating gradients at different reactor heights (red lines).Relevant experimental cells are open blue squares and arrows indicates species exchange between cells.The specific region of the substrate is green with its own cells interacting with the plasma phase.The dashed line is displaying the symmetry axis in this case and (B) to the corresponding full simulation box with communicating primary simulation boxes.Arrows gives the possible communications between the cells for covering all the experimental situations.Only a few are present for clarity.

FIGURE 2 |
FIGURE 2 | Schematics of the collision between two species (blue discs) both in simulations and in experiments and definitions of the interaction distance ρ, the distances d exp , d sim and times t exp , t sim between two collisions, respectively for experiments and simulations.v is the velocity which should be the same in experiments and simulations.

FIGURE 3 |
FIGURE 3 | Schematics of the simulation cell with the initial positions of all species involved in sputtering plasma simulation.

FIGURE 4 |
FIGURE 4 | Snapshot views of (a-c) the present simulation after 2 ns calculation time: with the substrate 12 × 12 × 6 nm 3 (a) without the substrate top (b) and side (c) views.Gas phase Ar close to the substrate are removed for clarity.And (d-f) from calculations presented in Xie et al. [44], on a 6 × 6 × 6 nm 3 similar porous substrate (d) without the substrate top (e) and side (f) views.The scale bar on the top of the figures holds for all pictures.The number of Pt atoms per unit area is the same.Blue and brown colored atoms are C and Pt respectively.

TABLE 1 |
Different plasma factors that can be accounted for in MD simulations.