Movaffaq Kateb

Jon Tomas Gudmundsson

Pascal Brault

Andrei Manolescu

Snorri Ingvarsson

On the role of ion potential energy in low energy HiPIMS deposition: An atomistic simulation

We study the effect of the so-called ion potential or non-kinetic energies of bombarding ions during ionized physical vapor deposition of Cu using molecular dynamics simulations. In particular we focus on low energy high power impulse magnetron sputtering (HiPIMS) deposition, in which the potential energy of ions can be comparable to their kinetic energy. The ion potential, as a shortranged repulsive force between the ions of the film-forming material and the surface atoms (substrate and later deposited film), is defined by the Ziegler-Biersack-Littmark potential. Analyzing the final structure indicates that, including the ion potential leads to a slightly lower interface mixing and fewer point defects (such as vacancies and interstitials), but resputtering and twinning have increased slightly. However, by including the ion potential the collision pattern changes. We also observed temporary formation of a ripple/pore with 5 nm height when the ion potential is included. The latter effect can explain the pores that have been observed experimentally in HiPIMS deposited Cu thin films by atomic force microscopy.

Introduction

Over the past few decades, high power impulse magnetron sputtering (HiPIMS) has attracted significant attention among the variety of ionized physical vapor deposition (PVD) techniques [START_REF] Helmersson | Gudmundsson, Ionized physical vapor deposition (IPVD): A review of technology and applications[END_REF][START_REF] Gudmundsson | The high power impulse magnetron sputtering discharge[END_REF]. In dc magnetron sputtering (dcMS) low plasma density (10 15 -10 17 m -3 [START_REF] Gudmundsson | Physics and technology of magnetron sputtering discharges[END_REF]) and low ionization fraction in the deposition flux (<10% [START_REF] Christou | Ionization of sputtered material in a planar magnetron discharge[END_REF]) is dictated by the thermal load at high powers, which might be destructive to the cathode target. As a result, the majority of the film forming species arriving at the substrate surface in dcMS are electrically neutral and the majority of the ions arriving at the substrate are ions of the working gas. In HiPIMS the overheating issue of the cathode target is resolved by applying unipolar high power pulses with a low duty cycle and at low repetition frequency, while maintaining similar time averaged power as for dcMS [START_REF] Helmersson | Gudmundsson, Ionized physical vapor deposition (IPVD): A review of technology and applications[END_REF][START_REF] Gudmundsson | The high power impulse magnetron sputtering discharge[END_REF]. Therefore a peak electron density that is 2 -3 orders of magnitude higher than in dcMS, is achieved in the vicinity of the cathode target [START_REF] Gudmundsson | Spatial and temporal behavior of the plasma parameters in a pulsed magnetron discharge[END_REF][START_REF] Bohlmark | Spatial electron density distribution in a high-power pulsed magnetron discharge[END_REF][START_REF] Meier | First measurements of the temporal evolution of the plasma density in HiPIMS discharges using THz time domain spectroscopy[END_REF]. Due to the high electron density the ionization mean free path of the sputtered species becomes shorter and the ionization probability increases and a significant fraction of the ions reaching the substrate are ions of the film-forming species. As a result HiPIMS deposition presents smoother [START_REF] Magnus | Morphology of TiN thin films grown on SiO 2 by reactive high power impulse magnetron sputtering[END_REF][START_REF] Sarakinos | Process characteristics and film properties upon growth of TiO x films by high power pulsed magnetron sputtering[END_REF][START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF], denser [START_REF] Samuelsson | On the film density using high power impulse magnetron sputtering[END_REF][START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF], and less defective [START_REF] Alami | Ion-assisted physical vapor deposition for enhanced film deposition on non-flat surfaces[END_REF][START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF] coatings, compared to dcMS deposition. Further, by applying a substrate bias, the bombarding energy of the ions of the film forming material can be tuned in order to achieve desired film properties such as the film texture and grain size [START_REF] Hajihoseini | Effect of substrate bias on properties of HiPIMS deposited vanadium nitride films[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF].

Molecular dynamics (MD) simulations have been shown to be a promising tool in terms of revealing atomistic mechanisms that contribute to various PVD techniques including dcMS and HiPIMS [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. It is worth mentioning that MD simulations of sputter deposition are divided into two categories: (i) modeling the vicinity of the substrate, the so called deposition side, and (ii) modeling the energetic working gas ions bombarding a surface or the target side. So far the target side simulations are focused on the sputter yield and energy distribution of ejected atoms from the target [START_REF] Brault | Molecular dynamics simulations of platinum plasma sputtering: A comparative case study[END_REF][START_REF] Kammara | Reconsideration of metal surface sputtering due to bombardment of high-energy argon ion particles: a molecular dynamics study[END_REF]. Recently, Brault [START_REF] Brault | Multiscale molecular dynamics simulation of plasma processing: Application to plasma sputtering[END_REF] included the working gas between target and deposition sides within a multi-scale MD simulation. He suggested that the effect of the gas phase in sputter deposition is limited to the narrowing of the energy distribution as atoms travel towards the substrate. It has also been shown that one can model the kinetic energy of adatoms at desired pressure without considering a working gas in the simulation [START_REF] Xie | Molecular dynamics simulations of clusters and thin film growth in the context of plasma sputtering deposition[END_REF]. Thus, it is reasonable to only model the deposition side using a proper energy distribution within the deposition flux.

Regarding HiPIMS deposition, MD simulations have indicated that both ionization fraction and substrate bias change the microstructure and film-substrate intermixing [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. Hence, the variation of the film texture by bias voltage that has been observed experimentally [START_REF] Cemin | Benefits of energetic ion bombardment for tailoring stress and microstructural evolution during growth of Cu thin films[END_REF][START_REF] Nedfors | Effect of synchronized bias in the deposition of TiB 2 thin films using high power impulse magnetron sputtering[END_REF][START_REF] Hajihoseini | Effect of substrate bias on properties of HiPIMS deposited vanadium nitride films[END_REF] can be attributed to the change of wetting due to interface mixing. An alternative explanation is based on thermal spikes [START_REF] Johnson | When is thermodynamics relevant to ion-induced atomic rearrangements in metals?[END_REF] or the atomic scale heating [START_REF] Musil | Hard and superhard nanocomposite coatings[END_REF] models of impinging ions. However, an MD simulation in the absence of ions is inconsistent with these models, as pointed out by Gilmore and Sprague [START_REF] Gilmore | Molecular-dynamics simulation of the energetic deposition of Ag thin films[END_REF]. Even energetic deposition (10 -40 eV) failed to make connection between MD results and the above mentioned models [START_REF] Sprague | Molecular dynamics simulations of filmsubstrate interface mixing in the energetic deposition of fcc metals[END_REF]. On the other hand, Müller [START_REF] Müller | Monte Carlo calculation for structural modifications in ionassisted thin film deposition due to thermal spikes[END_REF] showed that low energy ion impacts (ions of the noble working gas) can generate thermal spikes that lead to structure modification. Later, Anders [START_REF] Anders | Atomic scale heating in cathodic arc plasma deposition[END_REF] explained that the non-kinetic energy of ions, which will be called ion potential here after, cannot be ignored, and introduced several energy terms of an ion:

E i = E 0 i + ne∆V + E c + E iq + E exc + n Q=1 E Q , (1)
where E 0 i is the initial kinetic energy of the ion of the sputtered species, n is the ion valance, ∆V is the potential difference between the substrate and the plasma, E c is the cohesive energy, E iq is the image charge interaction energy, and E exc is a contribution that only applies to ions with excited electron(s). Finally, E Q is a cumulative ionization energy i.e. the summation of the first, second, ... and up to the n-th ionization energy to be considered for multiply charged ions. Note that the terms E 0 i and E c apply both to neutral atoms and ions. The term ne∆V takes into account the possibility to control an ionized flux using a substrate bias. Terms E c , E exc , E Q define the ion potential indicated by Anders [START_REF] Anders | Atomic scale heating in cathodic arc plasma deposition[END_REF][START_REF] Anders | Film deposition by energetic condensation[END_REF], which must be treated carefully in order to achieve a realistic simulation of an ionized PVD.

In an MD simulation E 0 i can be defined in terms of initial velocity. Since the simulated vacuum near the substrate is smaller than the substrate sheath (< 1 mm), it is not possible to model acceleration/deceleration of ions due to the substrate bias. Instead, it is convenient to consider the initial velocity corresponding to the summation of E 0 i and ne∆V [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. Since we are modelling the vicinity of the substrate, the above mentioned initial velocity refers to the velocity prior to impact.

The term E c is normally associated with the interatomic potential i.e. an individual atom in the vacuum has zero potential while the potential of an atom in the bulk is equal to E c . The image charge term is caused by the attraction between ions and image charges (electrons) that appears within a few nanometer depth from the surface [START_REF] Anders | Film deposition by energetic condensation[END_REF]. According to the over-barrier model [START_REF] Burgdörfer | Above-surface neutralization of highly charged ions: The classical over-the-barrier model[END_REF], for a low valence ion that approaches a metallic surface, neutralization occurs at a critical distance from the surface. This allows ignoring the E iq during deposition on a metallic surface [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. Both E exc and E Q can be taken into account using the electron force field (EFF) [START_REF] Su | Excited electron dynamics modeling of warm dense matter[END_REF] which is available for s and p electrons [START_REF] Jaramillo-Botero | Large-scale, longterm nonadiabatic electron molecular dynamics for describing material properties and phenomena in extreme environments[END_REF]. However, utilizing EFF is computationally intensive and increases the system size for taking into account electrons.

An interesting possibility is modeling the ion potential as an extra short-ranged repulsive interatomic potential [START_REF] Byggmästar | Effects of the short-range repulsive potential on cascade damage in iron[END_REF][START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. The term extra means that we already included an interatomic potential due to E c that might accurately describe forces around the equilibrium interatomic distances but not short-ranged ones [START_REF] Byggmästar | Effects of the short-range repulsive potential on cascade damage in iron[END_REF]. This method was first utilized by Müller [START_REF] Müller | Monte Carlo calculation for structural modifications in ionassisted thin film deposition due to thermal spikes[END_REF] to model working gas ions. He showed that bombarding a growing film by ions of the working gas increases film density and smoothness [START_REF] Müller | Ion-beam-induced epitaxial vapor-phase growth: A molecular-dynamics study[END_REF][START_REF] Müller | Stress and microstructure of sputter-deposited thin films: Molecular dynamics investigations[END_REF]. Also, a higher ratio of working gas ions to neutral atoms was found to reduce the number of voids [START_REF] Müller | Monte Carlo calculation for structural modifications in ionassisted thin film deposition due to thermal spikes[END_REF]. Later, Fang et al. [START_REF] Fang | Effect of gas impurity and ion bombardment on stresses in sputter-deposited thin films: A molecular-dynamics approach[END_REF] applied this model to the film forming ions and obtained a negligible effect from the working gas ions. Recently, Byggmästar et al. [START_REF] Byggmästar | Effects of the short-range repulsive potential on cascade damage in iron[END_REF] applied the model for a realistic simulation of cascade damage. They showed that the Ziegler-Biersack-Littmark (ZBL) potential [START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF]Chap. 2] can compensate the underestimation of the repulsive range by the embedded atom method (EAM) [START_REF] Daw | Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals[END_REF][START_REF] Daw | Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[END_REF]. As a result, the threshold displacement energy and many-body repulsion at short interatomic distances can be precisely modeled.

More recently, we utilized a combination of the ZBL potential and the EAM for a simulation of deposition with partially and fully ionized flux of sputtered species [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF]. It has been shown experimentally that surface roughness decreases as higher ionization fraction is utilized [START_REF] Magnus | Morphology of TiN thin films grown on SiO 2 by reactive high power impulse magnetron sputtering[END_REF][START_REF] Sarakinos | Process characteristics and film properties upon growth of TiO x films by high power pulsed magnetron sputtering[END_REF]. This occurs through the so-called bi-collision event, that causes a localized amorphization followed by recrystallization. The latter was not observed in the previous MD simulation studies without considering ions [START_REF] Houska | Quantitative investigation of the role of high-energy particles in Al 2 O 3 thin film growth: A molecular-dynamics study[END_REF][START_REF] Xie | Molecular dynamics simulations of clusters and thin film growth in the context of plasma sputtering deposition[END_REF]. We could also observe in the simulation results the formation of twinning during epitaxial growth of Cu using HiPIMS [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF], that had been reported in previous experiments [START_REF] Cemin | Epitaxial growth of Cu (001) thin films onto Si (001) using a single-step HiPIMS process[END_REF][START_REF] Kateb | Effect of atomic ordering on the magnetic anisotropy of single crystal Ni 80 Fe 20[END_REF]. Then, we showed that tuning the kinetic energy of ions through a substrate bias allows minimizing the creation of defects [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. A similar trend of polycrystalline to epitaxial transition using a substrate bias has also been observed experimentally [START_REF] Cemin | Benefits of energetic ion bombardment for tailoring stress and microstructural evolution during growth of Cu thin films[END_REF]. Moreover, it has been shown that energetic deposition (e.g. in highly biased HiPIMS or cathodic arc deposition) causes re-sputtering of deposited film and substrate, a process known as potential sputtering [START_REF] Winter | Interaction of slow HCI with solid surfaces: What do we know, what should we know ?[END_REF]. This effect explains the reduction of Cu deposition rate with increased substrate bias observed experimentally [START_REF] Cemin | Benefits of energetic ion bombardment for tailoring stress and microstructural evolution during growth of Cu thin films[END_REF].

The above mentioned examples clearly indicate the usefulness and reliability of describing the ion potentials using the ZBL term. The aim of the present manuscript is to study the effect of the ionic potential itself during the bombardment of a solid. We compare two cases, with and without considering the ZBL potential between ion and surface (substrate and later deposited film), and discuss changes in the film properties. We consider copper as the material of choice and assume a fully ionized flux of sputtered species to represent the HiPIMS process. This assumption is based on experimental findings from HiPIMS discharges with Cu target. For a Cu target, during HiP-IMS operation the ionization fraction of the sputtered flux up to 70% has been reported [START_REF] Kouznetsov | A novel pulsed magnetron sputter technique utilizing very high target power densities[END_REF]. We can also ignore the working gas ions because: (i) using HiPIMS, the working gas and metal ions arrive at the substrate at different times [START_REF] Macák | Ionized sputter deposition using an extremely high plasma density pulsed magnetron discharge[END_REF] and (ii) in the case of Cu it has been reported that up to 92% of ions arriving at the substrate are ions of the target material [START_REF] Vlček | Ion flux characteristics in high-power pulsed magnetron sputtering discharges[END_REF]. We neglect the contribution of E exc , as there is a low probability of finding excited ions far from the target [START_REF] Anders | Film deposition by energetic condensation[END_REF]. For Cu the first ionization energy is also relatively low (7.73 eV), and it can also be neglected compared to the ion kinetic energies [START_REF] Anders | Atomic scale heating in cathodic arc plasma deposition[END_REF].

Method

MD simulations [START_REF] Allen | Computer Simulation of Liquids[END_REF] were performed using the LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF][START_REF] Plimpton | Computational aspects of many-body potentials[END_REF] package 1 . The EAM force field was employed to model the interactions of film/substrate atoms. In the EAM, the total potential energy of atom i (E i ) is described by [START_REF] Daw | Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals[END_REF][START_REF] Daw | Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[END_REF]

E i = F i (ρ i ) + 1 2 i j φ i j (r i j ) , (2)
where F i is the embedding energy of atom i into electron density ρ i and φ i j is a pair potential interaction of atom i and j at distance r i j . Note that the electron density ρ i itself depends on electron density of neighboring atoms ρ i j so that ρ i = j ρ i j (r i j ).

The EAM potential relies on the same principles as the density functional theory and has been successfully applied to determine equilibrium and near-equilibrium behaviour of various metals and alloys. However, energetic deposition might push the system far away from its equilibrium. To describe the shortranged forces acting when an energetic atom travels inside the lattice we utilized the ZBL potential that has already been implemented in many force fields such as EAM [START_REF] Kammara | Reconsideration of metal surface sputtering due to bombardment of high-energy argon ion particles: a molecular dynamics study[END_REF][START_REF] Byggmästar | Effects of the short-range repulsive potential on cascade damage in iron[END_REF], Tersoff [START_REF] Bellido | Molecular dynamics simulations of ionbombarded graphene[END_REF] and deep-learning [START_REF] Wang | Deep learning inter-atomic potential model for accurate irradiation damage simulations[END_REF]

V(r i j ) = Z i Z j e 2 4πε 0 r i j Φ r i j a , (3)
where the Z i and Z j are the atomic numbers of species i and j that belong to the Coulomb term and e and ε 0 stand for elementary charge and vacuum permitivity, respectively. Note that the Coulomb term includes the atomic numbers in order to take into account very short interatomic distances (<2 Å) where electronelectron repulsion is dominating. Various ion potentials, such as Neilson, Brinkman and Firsove [50, Chap. 2], consider a similar Coulomb term while their screening function (Φ) might change.

The term Φ was introduced by Bohr [START_REF] Bohr | Atomic interaction in penetration phenomena[END_REF] and later modified by several authors. Finally, we apply the universal screening function Φ of Ziegler et al. [START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF] as discussed in our earlier work [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. The interested reader is referred to Townsend et al. [START_REF] Townsend | Ion implantation, sputtering and their applications[END_REF]Chap. 2] for comparison of various ion potential descriptions for the case of Cu. We considered two scenarios of ion-surface (substrate, and later deposited film) interaction: (i) EAM/ZBL and (ii) only EAM i.e. with and without ion potential, respectively. The repulsive nature of the ZBL potential may also result in exaggerated etching at high deposition rate. This issue can be solved by a checking mechanism, i.e. if an ion stays on the surface or implants into sublayers and remains there for 1 ps, it is considered to be an atom belonging to the film and its ZBL potential with respect to film atoms will be turned off. This 1 ps is consistent with the time scale for the electrons and the lattice to reach an equilibrium state after the localized heating due to high energy ion impacts [START_REF] Anders | Atomic scale heating in cathodic arc plasma deposition[END_REF].

It is also worth mentioning that the ion-ion interaction in the flux was modeled via a hybrid approach based on EAM and ZBL potentials. This might seem unnecessary since the flux travel towards the substrate is conventionally assumed to be collision-less. The latter is only true within a fixed energy flux, otherwise the energy distribution allows gas phase collisions [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. Recently, Brault [START_REF] Brault | Multiscale molecular dynamics simulation of plasma processing: Application to plasma sputtering[END_REF] demonstrated the gas phase clustering using a multi-scale MD simulation, i.e. when the deposition flux travels a distance similar to that in a typical experiment, with gas phase collisions treated properly.

The substrate was assumed to be a single crystal Cu, consisting of 16 (111) planes with 77×90 Å 2 dimensions. This makes the 111 orientation to be the growth direction. The substrate was divided into 3-layers as proposed by Srivastava et al. [START_REF] Srivastava | Anisotropic spread of surface dimer openings in the initial stages of the epitaxial growth of si on si {100}[END_REF]. A fixed monolayer is set at the bottom to prevent the other layers from moving after collisions, then three monolayers are set as a thermostat layer to control the heat dissipation and prevent melting, and the rest of the substrate is the surface layer. The initial velocities of substrate atoms were defined randomly from a Gaussian distribution, to mimic a temperature of 300 K, and the substrate energy was minimized prior to relaxation.

For both cases studied, the deposition flux consisting of 22000 ions was introduced at a distance of 150 nm above the substrate surface. The initial velocities of the ions towards the substrate were assigned randomly, with a flat distribution corresponding to the 0.1 -40 eV energy range. The flat energy distribution may seem a rather crude approximation. The energy distribution within deposition flux has been measured experimentally for both dcMS and HiPIMS. The original measurement of the ion energy distribution for Ti + ions in HiPIMS operation indicates a flat distribution up to 40 eV which then levels off [START_REF] Bohlmark | The ion energy distributions and ion flux composition from a high power impulse magnetron sputtering discharge[END_REF], while for Cu + it levels off around 20 eV [START_REF] Vlček | Ion flux characteristics in high-power pulsed magnetron sputtering discharges[END_REF]. The process of introducing species was a single ion every 0.1 ps, which produces an equal deposition rate in both cases. One may think this will generate deposition rate several orders of magnitude larger than in a typical experiment. We would like to remark that, considering the pulse duration, the instantaneous deposition rate of HiPIMS is 1 -2 orders of magnitude higher than for dcMS [START_REF] Kateb | Comparison of magnetic and structural properties of permalloy grown by dc and high power impulse magnetron sputtering[END_REF]. Besides, experimental measurements indicate a variable intensity of the ion flux during the pulse, reaching its peak within a fraction of the pulse [START_REF] Greczynski | Paradigm shift in thin-film growth by magnetron sputtering: From gas-ion to metal-ion irradiation of the growing film[END_REF]. Therefore, for HiPIMS the instantaneous deposition rate is typically several orders of magnitude higher than that of e.g. dcMS. However, the time scale that is achieved in MD simulation is generally on the order of tens of ns. Thus, MD simulations cannot capture the entire HiPIMS pulse which is normally 50 -400 µs long [START_REF] Gudmundsson | The high power impulse magnetron sputtering discharge[END_REF].

The time integration of the equation of motion was performed using the Verlet algorithm [START_REF] Verlet | Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF][START_REF] Kateb | Comparison of fracture behavior of sharp with blunt crack tip in nanocrystalline materials by molecular dynamics simulation[END_REF] with a timestep of 5 fs and sampling from the microcanonical (NVE) ensemble. The Langevin thermostat [START_REF] Schneider | Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions[END_REF] was only applied to the thermostat layer with a damping of 5 ps for a total time of 2.5 ns. The damping defines the timescale for resetting the temperature and that generates a heat flow towards thermostat layer.

In order to study the film-substrate interface quantitatively, we used the partial radial distribution function, g i j (r) introduced by Ashcroft and Langreth [START_REF] Ashcroft | Structure of binary liquid mixtures. I[END_REF], that was originally proposed for homogeneous binary mixtures. Here we have adapted it to the case of interface mixing, with the assumptions explained elsewhere [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF].

We considered the common neighbor analysis (CNA) method [START_REF] Steinhardt | Bond-orientational order in liquids and glasses[END_REF][START_REF] Kelchner | Dislocation nucleation and defect structure during surface indentation[END_REF][START_REF] Tsuzuki | Structural characterization of deformed crystals by analysis of common atomic neighborhood[END_REF] to characterize the local structure. CNA allows us to distinguish between face centered cubic (fcc) and hexagonal close-packed (hcp) structures which is of practical importance in defect analysis [START_REF] Kateb | Size and shapedependent melting mechanism of Pd nanoparticles[END_REF][START_REF] Azadeh | Determining phase transition using potential energy distribution and surface energy of Pd nanoparticles[END_REF]. This is achieved by de-termining the angles between the nearest neighbors, and thus the fcc and hcp lattices, having central and mirror symmetry in (111) planes, respectively, can be distinguished by a slight angle difference. The Ovito package 2 was used to generate atomistic illustrations and its Python interface was used for post-processing of CNA and g i j (r) [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with ovito-the open visualization tool[END_REF]. It is worth mentioning that Ovito allows utilizing an adaptive CNA, demonstrated by Stukowski [START_REF] Stukowski | Structure identification methods for atomistic simulations of crystalline materials[END_REF][START_REF] Stukowski | Computational analysis methods in atomistic modeling of crystals[END_REF] with variable cutoff calculation. This means that it is not sensitive to the choice of cutoff for finding the nearest neighbors. However, even the adaptive CNA characterizes atoms located at the surface as disordered, due to the lack of symmetry.

Furthermore, we studied vacancy and interstitial defects using the Wigner-Seitz (WS) method [START_REF] Nordlund | Defect production in collision cascades in elemental semiconductors and fcc metals[END_REF], based on the comparison of a current structure with an initial or a perfect lattice as reference. The method generates a primitive cell volume for each atom in the reference lattice and counts current atoms inside each volume. Thus, empty and over-occupied volumes translate to vacancy and interstitial defects.

Results

Intermixing

Figure 1 shows the cross-section of the film/substrate layers, obtained with and without considering the ion potential in the simulation. In both cases, severe intermixing is obtained which is a characteristic of low biased HiPIMS deposition process [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. However, a relatively thicker transition region is obtained when the ion potential is included in the simulation. Table 1 summarizes the results quantitatively. In both cases almost a 40 Å thick film is deposited. The exact values are presented by Z max film in Table 1. There is a deviation between the number of deposited species (22000) and the number of film atoms, N film . There are also losses in the substrate atoms indicated by ∆N sub . These quantities indicate a slightly increased resputtering when considering the ion potential. This difference is referred to as potential sputtering effect [START_REF] Winter | Interaction of slow HCI with solid surfaces: What do we know, what should we know ?[END_REF]. The value of Z min film refers to depth of the lowermost film atom which is about 26 Å below the initial substrate surface for both cases. We also located the topmost substrate atom in the film at the position denoted by Z max sub , which is at about 4 nm distance from original substrate surface. Assuming Z max sub -Z min film represents the interface thickness, a value of ∼60 Å is obtained with a negligible difference i.e. it is 4 Å larger with ion potential included in the simulation.

Microstructure

The CNA analysis of films deposited with and without considering the ion potential is shown in figure 2. Various defects can be observed and characterized in both cases. The presence of a hcp layer in a fcc material is generally referred to as a stacking fault, or a twinning when it is large enough. Vacancies can be detected using CNA as empty sites enclosed by 12 disordered atoms. A tetrahedral or octahedral configuration of 2 Ovito website, http://ovito.org/, Version 3.0.0-dev794 2 indicates the original substrate surface that is used to distinguish between the film and substrate defects.

It can be seen that, in both cases, low energy HiPIMS deposition induces a considerable number of various defects into the substrate. Introducing defects into the substrate is the characteristic of ionized PVD and does not occur during e.g. evaporation when almost all the film-forming species are neutral [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Stencil growth of metallic nanorod: An atomistic simulation[END_REF]. However, as previously demonstrated [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF], high ionization fraction of HiPIMS allows us to tune the kinetic energy of the flux using the substrate bias and efficiently minimizing the defects. WS analysis also indicates a considerable number of point defects, both in the film and in the substrate. Note that, however, the density of vacancies is more uniform in the film and substrate, while interstitials are mostly located in the films. This means that both films here present a more compressive stress than their substrates.

Regarding the effect of ion potential on the microstructure, it can be seen that without the ion potential more point defects are introduced both in the film and in the substrates. However, the ratio of hcp atoms is higher when the ion potential is taken into account.

Temporal behaviour

Figure 3 illustrates g i j (r) as a function of the elapsed time during the simulated deposition with and without the ion potential. The color map indicates the intensity of the normalized g i j (r) i.e. the density of film-substrate bonds. The major peak indicates the distance between the first nearest neighbors (1NN). That is 2.49 Å at the early stage of deposition, and becomes 2.55 Å with the formation of the film. Such an increase in the 1NN generally indicates an increase in the coordination number, or, here, per atom film-substrate bonds. The disturbances in the g i j (r), which are more evident in the major peaks, mean the redistribution of atoms due to collisions. Each collision shifts the 1NN peak to lower distances, that indicates a compression of the system. Then, the 1NN peaks relax to their original position. The arrows point to four major collisions that cause a clear shift in the 1NN peak. For a flat interface it is expected that g i j (r) tends to zero as the film thickness increases. This is because in a flat interface the number of film-substrate bonds remains constant while the total number of bonds grows. In both cases here, the intensity of the major peak in general increases. This is an indication of a continuous interface mixing and the development of more film-substrate bonds.

Next, figure 4 shows the time evolution of the microstructure during deposition. The arrows here correspond to collisions indicated in figure 3. Looking into the fcc percentage, it can be seen that after the indicated collisions the system requires some time for recovery. The only exception is collision 1 in figure 4(a), that occurs at the early stage of deposition. The collisions indicated by arrows produce a severe amorphization followed by a step change in the hcp percentage. We observed that the rest of the collisions can only produce point defects. To our knowledge Klokholm [START_REF] Klokholm | Intrinsic stress in evaporated metal films[END_REF] was the first to suggest the existence of a disordered liquid-like material buried beneath the advancing film surface during evaporation. He believed that, during the evaporation, the energy of adatoms (0.1 -4 eV) is not high enough to cause the recovery of sub-surface layers. Instead, we used ions and much higher energies (0.1 -40 eV), that allow for amorphous to crystalline phase change. However, even this range of energies is not sufficient for the formation of equilibrium fcc structure to proceed without defects. The latter can be achieved only at higher energies, like 60 -100 eV [START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF].

We would like to recall that we used similar initial velocities in both cases studied. We found that when the ion potential is considered the collision pattern changes. Note that without the ion potential (see figure 4(b)) the collisions indicated by arrows are more regularly spaced compared to the case with the ion potential included in the simulation. This is mainly because the presence or the absence of the ionic potential leads to a variation of the kinetic energy in vicinity of the surface.

Figure 5 shows a temporary surface roughness that appears after collision 3 in figure 4 (a). This ripple/pore remains almost unchanged until collision 4. The formation of similar ripples, with significant surface roughness, has been observed after ion irradiation of Si surface under an angle using kinetic Monte Carlo (KMC) simulations within binary collision approximation [START_REF] Liedke | Surface morphology and interface chemistry under ion irradiation-simultaneous atomistic simulation of collisional and thermal kinetics[END_REF], hybrid MD/KMC [START_REF] Yang | Kinetic monte carlo simulation of selforganized pattern formation induced by ion beam sputtering using crater functions[END_REF], and also experimentally [START_REF] Madi | Multiple bifurcation types and the linear dynamics of ion sputtered surfaces[END_REF][START_REF] Keller | Simultaneous formation of two ripple modes on ion sputtered silicon[END_REF]. While these studies were focused on ions incident under an angle, e.g. 80 • , we considered the initial velocity of flux perpendicular to the substrate surface. However, the appearance of a pore indicates an elongation that is a characteristic of oblique incidence. Thus, the deflection of the energetic ions was achieved due to the repulsion between the ions and the surface. We did not observe the formation of any pore without considering the ion potential. This is due to the fact that considering ion potential would lead to a more realistic collision behaviour, such as possibility of deflection, while the average behaviour or overall behaviour for longer timescale may seem similar (cf. figures. 1-4). We would like to remark that the formation of open pores has been verified by atomic force microscopy of Cu thin films deposited by HiPIMS [START_REF] Cemin | Benefits of energetic ion bombardment for tailoring stress and microstructural evolution during growth of Cu thin films[END_REF]. To the best of our knowledge, in this work we have shown for the first time that such effects can occur due to the ionic potential. The formation of open pores at the surface has also been reported for relatively harder materials, such as CrN and NbN deposited by HiPIMS [START_REF] Biswas | Effect of substrate bias voltage on defect generation and their influence on corrosion and tribological properties of HiPIMS deposited CrN/NbN coatings[END_REF]. Thus, our findings are of general interest for variety of thin film deposition processes.

It is worth mentioning that earlier we reported enhanced surface roughness obtained by utilizing a flux of neutrals [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF][START_REF] Kateb | Effect of substrate bias on microstructure of epitaxial film grown by HiPIMS: An atomistic simulation[END_REF]. Such a surface roughening were caused by so-called steering effect [START_REF] Van Dijken | Steering-enhanced roughening during metal deposition at grazing incidence[END_REF]. The latter causes surface roughness by preferential landing of arriving neutrals on top of islands due to the long-ranged attractive forces. However, the ion potential here was applied through the ZBL potential which is of short-ranged repulsive nature. Besides, the pore shown in figure 5 of the present manuscript, appears immediately after collision 3 (figure 4(a)) and disappears after collision 4 while surface roughnesses observed using a flux of neutrals were formed gradually.

Discussion

A question that may rise here is when is it necessary to include ion potential? To answer the question we need to highlight couple of considerations. As mentioned earlier, most of the interatomic potentials, such as EAM are obtained by fitting some particular properties at equilibrium condition. Thus, low energy collisions, e.g. during deposition of neutral adatoms or partially ionized flux, that does not push the system out of equilibrium can be described without including ion potential. For high ionization flux fractions, the probability of collisions that may go far from equilibrium increases and one may consider including ion potential. This also depends on the ions kinetic energy i.e. very high kinetic energies deposit much more energy into the film which then makes the influence of the ion potential negligible. On the other hand, in the low kinetic energy regime, there is lower probability of collisions with non-equilibrium behaviour. However, there are conditions where ion potential must be included such as when modelling the target side.

It is also worth mentioning that, conventionally ionized PVD methods have been modeled with energetic neutrals rather than ions (cf. Ref. [START_REF] Kateb | Role of ionization fraction on the surface roughness, density, and interface mixing of the films deposited by thermal evaporation, dc magnetron sputtering, and HiPIMS: An atomistic simulation[END_REF] and references there in). Now another question arises: can energetic neutrals compensate for the ion potential? The answer is that the ion potential is effective in a very short interatomic distance that does not occur in all collisions. Thus, the effect of the ion potential is very localized. However, energetic neutrals bombard everywhere on the surface and modify the structure more homogeneously.

Summary

In summary, the effect of the ion potential energy on a film deposited by highly ionized flux of the film-forming material, such as during HiPIMS deposition is studied using MD simulations. We considered two cases, with and without considering the ZBL repulsive potential between ions and surface. We have shown that there is a slight difference between the two cases, in terms of surface roughness, interface mixing, point defects, the microstructure, and the time evolution. In addition, we have shown that the ion potential changes the collision pattern. In particular, we have obtained a temporal ripple/pore with 50 Å height difference, that can explain experimental results.

Figure 1 :

 1 Figure 1: Cross-section of the film/substrate layers (a) with and (b) without considering ion potential. The film atoms are illustrated by larger red atoms, and substrate atoms are depicted smaller yellow atoms, for clarity. Note that the substrate atoms are shifted in front of film atoms for a better illustration of the transition region. The scale on the right shows height from the initial substrate surface.

Figure 2 :

 2 Figure 2: Demonstration of the microstructure using CNA for deposition (a) with and (b) without ion potential. The yellow dashed line indicate initial substrate surface. The fcc atoms are shown as semitransparent for illustration of internal defects.

Figure 3 :Figure 4 :

 34 Figure 3: Variation of g i j (r) during the deposition (a) with and (b) without ion potential. The colormap indicates normalized density of film-substrate atom pairs. Arrows indicate considerable shift in 1NN peaks due to major collisions.

Figure 5 :

 5 Figure 5: Temporary surface roughness due to the ion potential obtained (a) before and (b) after collision 3 (figure 4(a)), and (c) before collision 4.

  potentials. The ZBL potential is a modification of screened Coulomb potential[START_REF] Ziegler | The Stopping and Range of Ions in Matter[END_REF] Chap. 2]

Table 1 :

 1 N film denotes number of film atoms and ∆N sub refers to sputtered substrate atoms, and Z indicate maximum/minimum position of the film/substrate atoms with respect to the initial substrate surface.

	Method	with	without
	N film (atom) 21646 21733 Z min -25.8 -25.9 film (Å) Z max 42.8 44.0 film (Å) 202 147 ∆N sub (atom) Z max sub (Å) 40.8 36.6

Table 2 :

 2 The results from the WS analysis of the final structure with N V and N I being the number of vacancies and interstitial atoms, respectively. represent an interstitial. Further details on the point defects analysis by WS are summarized in Table2. The dashed line in figure

	Method	with without
	N V sub (atom) N I sub (atom) N V film (atom) 1066 911 427 N I film (atom) 1225	959 448 1268 1288

LAMMPS website, http://lammps.sandia.gov/, distribution 14-April-2018

Acknowledgements

This work was supported by the University of Iceland Research Funds and the Icelandic Research Fund, grants nos. 195943, 196141, 130029 and 120002023.


A Fourmont
email: adrien.fourmont@u-bourgogne.fr

) O Politano

S Le Gallet
email: sophie.le-gallet@u-bourgogne.fr

C Desgranges
email: clara.desgranges@safrangroup.com

F Baras

Reactivity of Ni-Al nanocomposites prepared by mechanical activation: a molecular dynamics study

High energy ball milling of metallic powders leads to high reactivity in the milled mixture. The reaction is often faster and starts at a lower temperature.

However, the mechanisms responsible for this high reactivity are not yet completely understood. The aim of this study is to evaluate one of the possible activating factors of this heightened reactivity: nano-scale mixing of the reagents. Molecular Dynamics was used to analyze the role of an amorphous Ni-Al mixing layer, mimicking the powder microstructure after milling, between two Ni layers. The impact of temperature and stoichiometry was investigated in relation to the formation of the B2-NiAl intermetallic compound. At low temperatures, pre-mixing does not seem to slow down the diffusion of Ni atoms in an amorphous Al region. Homogeneous nucleation was observed in this peculiar milled microstructure. These two phenomena explain why the nano-scale mixing observed experimentally after high energy milling is indeed an activating factor in the reactivity of metallic systems such as Ni-Al.

I. INTRODUCTION

Milling several compounds in order to mix them is a very old and common technique. However, high energy milling is relatively new to the domain of metallurgy 1 .

Despite its widespread use today, the high energy ball milling of metallic powders is not completely understood, especially in terms of the evolution of the microstructure at atomic scale during the milling or even the description of ball movement 2 . In the area of Self-Propagating High Temperature Synthesis (SHS), high energy ball milling is a well-known elaboration process. First, two compounds are milled together and then the milled powders are heated at a low temperature (up to 900K). When the ignition temperature is reached, the reaction starts and self-propagates. It has been observed that milling two powders leads to high reactivity and solid-state reactions, which translate into a lower ignition temperature and a faster reaction [3][4][5] . In this situation, the amount of the liquid phase is reduced, and the reaction occurs at a moderate temperature as compared to other SHS processes such as SHS in nanolaminates. One of the more widely studied systems is the Ni-Al system which forms the intermetallic Ni 50 Al 50 compound with an ordered B2 structure. To explain this high reactivity, Rogachev et al.

studied milled powders of Ni+Al by means of High-Resolution Transmission Electron

Microscopy (HRTEM) and suggested 6 activating factors 6 :

(i) decreasing apparent activation energy by the magnitude of the energy stored in reactants during MA;

(ii) increasing atomic mobility due to high concentration of the grain boundaries and defects created by plastic deformation;

(iii) formation of non-equilibrium meta-stable phases;

(iv) increase in the contact (interfacial) surface area between reactants;

(v) formation of freshly formed reactive surfaces; and

(vi) nano-scale mixing of the reagents.

These factors are nearly impossible to thoroughly study experimentally because of the time and space scales. Despite the recent success of innovative experiments investigating milling regimes in planetary ball mills 7 , today's technologies do not make in-situ microscopic observations possible. Thus, microscopic simulations such as Molecular Dynamics (MD) could provide a useful tool to describe both milling and SHS processes. The milling mechanism has been modelized in several ways in the literature.

The sliding of metallic blocks and the subsequent friction leading to mixing has been studied by several research groups [8][9][10] . They have observed mixing at the atomic scale in an amorphous zone and have noted that the mixing resulted from ballistic diffusion rather than classic thermal diffusion. These teams also observed the formation of structural defects resulting from this friction. Chen et al. studied the sliding of blocks of Ni and Al. They have reported that the mixing layer is mainly formed in the softer Al material 11,12 . Odunuga et al. and Delogu et al. evaluated another kind of simulation using cyclic deformations that mimic milling 13,14 . They studied forced chemical mixing driven by plastic deformation and observed that the mixing resulted from dislocation gliding and that it was similar to the mixing of a substance subject to a turbulent flow.

The reactivity of these simulated systems is to analyze because two main, interdependent aspects have to be considered: the creation of structural defects and the induced chemical mixing. Adding to the complexity of analysis is the fact that the importance of one aspect could overtake that of the other due to issues of scale.

Microscopic simulation of the reactivity of a milled system is a challenging undertaking because of the co-occurrence of many factors. Preliminary attempts to examine the reactivity of activated systems have been carried out in specific cases. der particles and the subsequent reactivity [15][16][17][18] . Here, the rise in temperature is a consequence of the impact, contrary to the case of a classic SHS process. This is a one step process: shock-induced chemistry. The authors demonstrated that in this case, compres-sion is associated with the sharp rise in temperature and that the subsequent reaction is dominated by interdiffusion. They also showed that porosities play a major role by acting as hot spots during the deformation 17 . Witbeck et al. have evaluated the effects of lattice defects and grain boundaries 19,20 . This work corresponds to an exploration of activating factor (ii). The authors observed that point defects greatly increased the reaction speed in SHS propagation, whereas grain boundaries acted as heterogeneous nucleation sites in the formation of the intermetallic compound B2-NiAl which further inhibited Ni diffusion. This study was performed in the context of nanolaminates where temperature increased tremendously (up to 2200K); indeed, temperature is the key factor in the reaction. However, these findings are not necessarily relevant to the case of milled systems, in which temperature and the resulting diffusion are not likely to be the main and (v): the effect of surfaces 21,22 . They have reported that the misfit strain decreased the activation energy of the reaction. They also observed that voids and porosity led to a slower mechanism compared to the case in which the two elements were in direct contact. The Weihs group has indirectly studied (vi) in the case of temperature driven mixing [23][24][25] . They noted that a pre-mixing step hindered diffusion and that a nano-scale gradient suppressed the homogeneous nucleation of B2-NiAl. They also observed that B2-NiAl readily nucleates and grows at the interface between an amorphous pre-mixed region and an f cc-Ni block. 26,27 . They analyzed the crystallization resulting from decreasing temperature. They observed that nucleation occurred in two steps: segregation of the Ni and Al into an amorphous but ordered structure, and, then, crystallization into the B2 compound. They observed that B2-NiAl homogeneous nucleation began at 900 K. Examining the case of a bilayer configuration, Baras et al. have reported that the crystallization of B2-NiAl resulted from a liquid phase 28,29 .

A systematic study of the effect of mixing (vi) in the framework of milled systems is still lacking. Indeed, two main parameters significantly differ from those which apply to studies on laminate systems: first, the temperature range is significantly lower and secondly, the reaction mainly occurs in solid-state. In addition, the microstructure of such systems includes more defects in the form of large amorphous zones generated during the milling process. It is nearly impossible to take into account all the significant parameters to model a real-world milled Ni/Al reactive materials at once. In the present work, we focused on the effect of pre-mixing in an amorphous system and on the resulting crystallization of B2-NiAl at low temperatures (below the melting temperature of Al). Indeed, a careful examination on the current literature on the topic led us to believe that this viewing angle was of primary importance 30 . This study modelized a two-step process in which the reactants are first milled and then heated, as is usually the case in such experiments.

II. SIMULATION DETAILS

To describe the Ni-Al system, we employed the classic Molecular Dynamics (MD) approach using the Embedded-Atom Method (EAM) potential developed by Purja Pun and Mishin 31 . This interatomic potential has been fitted to reproduce a large set of physical quantities determined experimentally or by ab-initio calculations. It therefore accurately describes the main physical properties of the pure elements (i.e. Ni and Al) as well as selected intermetallic compounds, namely B2-NiAl and L1 2 -Ni 3 Al. Despite its good transferability, this potential necessarily deals with a simplified system containing a limited number of phases as compared to a real Ni-Al system. Table I reports the melting temperature, cohesive energies and lattice parameters computed using this potential 32 .

The MD simulations were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) open-source software 33 and the atomic visualizations were effected via the Open Visualization Tool (OVITO) 34 .

Following on from our previous work 28 , we considered a layered Ni-Al-Ni system with the following characteristics. The reference system is made up of one layer of Al of 35 atomic planes in between two Ni layers of 36 atomic planes each along the z-direction (Fig. 1a). Both pure metals are in f cc configuration and the interface orientation is (001). The typical size of the simulation box is L x = L y ∼ 14 nm and L z = 20 nm.

The system is submitted to periodic boundary conditions in all directions. The initial distance δ between the Ni-and Al-atomic planes at the interface is 2.5 Å. The system is made up of 305 214 atoms (n(Ni) = 231 276 and n(Al) = 73 938). The inner layer is defined as the layer of Al plus 3 atomic planes of Ni interfaces on each side. The width in the x and y directions of the inner layer is chosen to be smaller than that of the outer Ni layers (L x (inner) < Lx(outer) and L y (inner) < L y (outer)). The thickness of the inner layer is 8.3 nm. This specific geometry including a narrower inner layer allows us to avoid artificial constraints due to the application of periodicity on both inner and outer layers during the construction of the system and its subsequent evolution. Each system is then denoted as 'REFxxK' where xx K refers to the temperature of the simulation (see Table II). II). The system contains 305 013 atoms with n(ML) = 85 220 in the mixing layer. For example, the mixing layer in the ML20% system contains 17 036 Ni atoms. In order to ensure that most of the mixing layer was truly amorphous at 300 K, we chose to consider the pair correlation function g(r) of the mixing layer as our main criterion. The g(r) functions of the inner layers of our different systems are displayed in Fig. 2. The pair correlation function of the mixing layer created in this way is similar to that of a liq-Al system (i.e.

melted Al) generated by heating f cc-Al at 2000 K over a long period: 3 ns. Hence, the amorphous character of the mixing layer is indeed mimicked by means of this rapid heating-cooling treatment. Note that the substitution ratio does not modify the degree of amorphisation. The configuration of the initial system with a mixing layer is depicted in Fig. 1b. In order to study the reactivity of such systems, both REF and MLXX% systems were studied at various temperatures in the isothermal-isobaric ensemble (NPT) in the range of 600 K-1200 K around the bulk melting temperature of aluminum.

Several indicators allow us to follow the evolution of the system at the microscopic level. If we calculate the local lattice structure according to the common neighbor analysis, we obtain the local configuration around a given atom. It is possible to distinguish between different configurations: bcc (1), f cc (2), hcp, and unknown (0). Each atom is thus labeled by the color of the species and by its local lattice structure. The potential energy per atom, which is very sensitive to the local environment of a given atom, was also evaluated. In addition, global indicators were followed during the evolution of the system: the atomic percentage of Ni atoms in the inner or mixing layer; the amount of bcc-atoms in at.% out of the maximum that could be formed

Ξ = 100 × n(bcc) 2 × n(Al) (1)
where n(bcc) and n(Al) are the numbers of B2and Al-atoms in the mixing layer (B2 atoms form a bcc lattice if we do not make the distinction between atomic species). If all Al-atoms are bcc, Ξ is close to 100 % and the phase is pure B2-NiAl. Ξ makes the comparison between several stoichiometries possible. A cluster analysis was developed in order to analyze the nucleation and growth of the intermetallic phase during evolution of the system. A cluster of atoms is defined as a set of neighboring bcc atoms. If the cluster is made up of at least 20 atoms, the cluster is counted as a nucleus.

III. RESULTS AND DISCUSSIONS

In a bilayered system Ni-Al-Ni, the typical evolution at an elevated temperature consists of several stages: the melting of aluminum, the progressive dissolution of Ni-atoms into the liquid layer and, eventually, the formation of B2-NiAl at interfaces. The corresponding scheme of reaction reads as follows:

Ni(s) + Al(l) → Ni(s) + AlNi x (l) → Ni(s) + AlNi x (l) + NiAl(s) (2)
This mechanism is named mosaic dissolution-precipitation 35 Next, we considered a system with a premixed layer that mimics an activated Ni/Al system after milling. The mixing layer contains 50 at.% of Ni and the bilayer is heated at 1200 K. The behavior of the ML50% system is actually very different from that of the REF1200K. There is no delay in the formation of B2-NiAl and the ratio Ξ reaches a higher value close to 50%. The decrease in unk-atoms (amorphous) is concomitant with the formation of B2-NiAl. As amorphous Al-and Ni-atoms are present in the mixing layer from the start, the reorganization of atoms takes place from the very beginning.

No prior melting was required. Note that the number of f cc-atoms increased slightly at the start of the phase transformation before reaching a plateau. We observed that the number of f cc-atoms is higher in the ML50% system as depicted in Fig. 9.

The main differences observed are certainly due to the fact that the melting and dissolution steps have already been performed in ML50%. However, it has been noted that a mixing layer hinders diffusion 23 . Indeed, according to Fick's Law, the concentration gradient decreases and diffusion slows down. In addition, diffusion of Ni atoms in a

Nickel amorphous zone is known to be faster than in a Nickel crystalline zone, as it was reported by Manukyan et al. for instance 36 ; in this case, the use of amorphous Ni as reactant was proved to be more efficient than crystalline Ni.

In order to evaluate the effect of melting/dissolution as compared to diffusion in B2formation, we studied the diffusion of Ni-atoms in both systems. Figure 4a depicts the evolution of the Ni concentration (at.%) in the mixing layer and the evolution of the B2 phase in the ML50% and REF systems. Ni concentrations were calculated in the whole mixing layer (solid lines) and in the 1 nm slice just above or below the interface (dashed lines) in order to evaluate the dissolution of Ni coming from the two f cc-Ni blocks. We also considered the ML20% system at 1200 K, where the inner layer is not saturated in Ni. The evolution of Ni-concentration can be correlated with the appearance of B2-nuclei as shown in Fig. 4b.

In ML50%, the Ni concentration does not evolve much since the inner layer is close to saturation in Ni 37 . There is nevertheless some limited diffusion at the interface. The only phenomenon which occurs is the ordering in B2. Nucleation proceeds very quickly until 0.5 ns at which point the seeds coalesce. In the end, there is one main grain of 34000 atoms. In a perfect B2 structure, the grain corresponds to a total volume of 410 nm 3 . Assuming the grain is a sphere, it has a radius of 4.6 nm.

In ML20%, the Ni concentration in the 1nm slice close to the interface undergoes a sharp increase until 0.5 ns, then continues to increase. The Ni concentration in the entire mixing layer increases more slowly and reaches a plateau value of approximately 35% at 3.5 ns. The inflexion point in the evolution of Ξ in ML20% corresponds to the plateau in Ni concentration. The formation of B2-NiAl seems to begin efficiently when the Ni concentration in the whole layer reaches a steady value. From 0.5 to 2 ns, Ni-atoms migrated from interfaces to the 'center' of the mixing layer and B2 seeds appeared.

After 1.5 ns, nucleation takes place with an increasing number of nuclei (Fig. 4b). After 3 ns, the number of nuclei decreased but some grains grew towards the 'center' of the mixing layer. At 8 ns, a limited number of rounded grains clung to the interfaces (Fig. 4c). In ML20%, phase formation proceeds through heterogeneous nucleation: no grains formed spontaneously in the bulk of the mixing layer.

The behavior of the REF sample is similar to that observed in ML20%, except that the Ni concentration reaches its steady value slightly later and this value is lower. Nucleation of B2-NiAl grains is shifted to 2 ns. In addition, fewer but larger grains are formed (biggest nucleus r = 3.2 nm). Small grains merged with larger ones, as there is some coalescence (see Supplementary Materials).

At high temperature (1200 K), the REF system and the ML system with a low Ni concentration exhibit similar behavior: dissolution of Ni close to the interface, heterogeneous nucleation, further dissolution of Ni to the 'center', growth. But in the case of the premixed layer, the formation of B2 occurs sooner and the system is more reactive. Higher Ni concentrations in the mixing layer lead to a shift in the start of nucleation. In addition, the final amount of B2-NiAl is larger in higher concentrations. This observation holds until 50%. For concentrations of Ni equal to 60% and higher, the f cc solid solution forms more rapidly at the interface and slows down the formation of B2-NiAl.

The Al cohesive energy/atom E 0 (Al) depends on the structure to which it belongs: the decrease in E 0 (Al) is certainly not the only factor responsible for the high reactivity observed in premixed systems. As evidence of this, the ML60% system does not form more B2 than the others.

As stated above, the purpose of this work was to evaluate the formation of B2-NiAl at low temperatures to mimic the ignition of milled Ni-Al which usually occurs experimentally at around 750 K 6 . The temperatures chosen for our simulations are higher in order to obtain results within an acceptable computing time. First, the actual Al melting temperature in the particular geometry was roughly estimated. Figure . 7a gives the number of unk-Al as a function of time in the reference system at different temperatures. At 800 K, most Al-atoms remain f cc. At 900 K, after 4 ns, all the aluminum has melted. In this geometry, the effective Al melting temperature seems to be in the range of [800 K -900 K], much lower than the bulk melting temperature of 1055 K (see Table I). However, simply meeting the criterion of a melted inner layer is not enough to ensure the formation of B2-NiAl. Indeed, the REF systems at 900 K up to 1100 K do not exhibit the formation of B2-NiAl (see Fig. 7b).

The evolution of the Ξ factor is plotted in Fig. 7c at different operating temperatures in the case of ML50%. In this pre-mixed system, B2-NiAl is formed even at 800 K, below the effective melting temperature. The (Al+Ni) solid is in an amorphous or supercooled state. Hence, for a high Ni concentration in the mixing layer, this configuration already exhibits solid state reaction at 800 K. Formation of B2-NiAl becomes more significant at 900 K. It is important to compare the mechanism of B2-NiAl formation at low (800 K and 900 K) and high temperatures (1200 K) at the different rates of premixing from 20% up to 50%. tion, the intermetallic B2-NiAl phase appears plane by plane. The atomic arrangement is characteristic of the (110) orientation in NiAl. In each plane, only four orientations in (110) planes were observed, due to the specific orientation relationships between bccand f cc-phases 29 . Thus, small disoriented grains separated by grain boundaries were observed. The situation is radically different in the ML50% system at 900K and 800K.

Seeds appear in the middle of the layer, far from interfaces. This is the signature of a homogeneous nucleation. In addition, the typical profile depicted in Fig. 8 shows that the inner layer remains compact and concave in a solid amorphous state. Wetting of the Ni substrate is observed, but its impact remains limited since the layer is already amorphous. This means that the homogeneous nucleation occurs as a solid reaction.

The behavior is similar at 800 K and 900 K, except for the fact that B2 is formed more rapidly at 900 K. These findings indicate the possibility of B2-NiAl formation at temperatures below the melting point of aluminum and could explain the enhanced reactivity of milled systems.

Experimentally, the intermetallic L1 2 -Ni 3 Al phase is sometimes observed in this kind of study. In our simulation work, it does not seem to be the case. The observation of homogeneous nucleation in a mixing layer at 900K leads to other questions. Is homogeneous nucleation caused solely by the 50% pre-mixed layer?

Would another premixing stoichiometry also lead to homogeneous nucleation? These questions are of great interest since an experimental milled system could be modelized as the sum of several zones, each with a different stoichiometry. Three compositions prone to homogeneous nucleation at 900 K were more closely studied: ML30%, ML40% and ML50%. ML50% at 900K exhibits a different behavior. This system shows a slower but steady nucleation rate. Grain size does not increase much, as revealed in the global evolution.

At the end, there are approximately 40 nuclei and the largest grain does not exceed the size of 2000 atoms (r = 1.8 nm).

In order to investigate the differences between behavior in ML50% and ML40%, representative grains in each system were followed over the course of their formation.

For selected systems, a significant time was chosen. This time corresponds to the complete formation of nuclei of reasonable size. From that time, all the BCC atoms in a selected nucleus are then selected. This observation is represented in Fig. 11b illustrates the rearrangement of atoms during homogeneous nucleation in HON40:

randomly distributed atoms diffuse and explore a large domain for 1.5ns before the rapid crystallization of the B2-NiAl phase. The situation is radically different in HON50, as shown in Figure . 11c. From the start, a fifth of the atoms form a seed of bcc-atoms that barely evolves up to 1.2ns. At that time, some new atoms adhere to the seed and then the rest of atoms abruptly become bcc. Crystallization begins sooner in the ML50% system, but since mobility is hindered, the atoms truly rearrange into a crystalline structure only after 1ns.

Fig. 12 depicts the evolution of the diffusion coefficient of Ni at 900K according to concentration. The maximum is reached at around 30% and the diffusion at 40% is 1.3 times higher than at 50%, according to these calculations. There are indeed different mechanisms at work in the ML40% and ML50% systems. The higher mobility in ML40% is due to the lower Ni concentration. The high reactivity of activated systems results from a compromise between diffusion and previously effected mixing.

IV. CONCLUSIONS

The bilayer configuration with a pre-mixed zone is indeed characterized by high reactivity. Here, reactivity was evaluated by measuring the amount of B2-NiAl formed during homogeneous heating. In this pre-mixed configuration, two steps usually performed over time by temperature effect are already achieved in the initial configuration.

The melting and amorphization of Al and the diffusion of Ni are already accomplished.

Conversely, diffusion is slower because of the pre-mixing. Taking into account the two conflicting phenomena, we nonetheless note that the pre-mixed configuration results in higher reactivity compared to the reference bilayer configuration. At the microscopic scale, it has been observed that pre-mixing decreases the cohesive energy / atom by as much as 40% for the pre-mixing with 50% Nickel.

At low temperatures below the melting temperature of Al, homogeneous and heterogeneous nucleation occur in concentrated pre-mixing zones. The presence of both kinds of nucleation leads to faster formation of higher amounts of B2-NiAl, whereas in systems where only heterogeneous nucleation occurs, the formation is slower. It is interesting to note that the equimolar stoichiometry 50% Al-50% Ni is not associated with the best reactivity at low temperature. The mobility of atoms indeed seems to be the predominant factor in the onset of reaction. A slightly less concentrated system (ML40%) is a better compromise to ensure higher reactivity.

In addition, in cases in which the system is not pre-mixed, the Ni/Al interfaces are rapidly saturated and the formation of B2 slows down the interdiffusion of Ni. The intermetallics formed at the interfaces seem to channel the Ni diffusion through grain boundaries. Thus, the formation at the center of the inner layer is possible but inhibited.

In conclusion, our simulations demonstrate that the nano-scale mixing of the reagents (factor (vi) as suggested by Rogachev et al.) is indeed an activating factor. Experimentally, the two main conclusions drawn in this work could still stand. At low temperature, pre-mixing does not seem to slow down the diffusion of the Ni atoms in an amorphous Al region. In addition, the homogeneous nucleation resulting from the peculiar milled microstructure would accelerate the general reaction of B2-NiAl formation. Interestingly enough, similar results were obtained in the case of irradiated films where a Ni/Al amorphous interlayer was generated and found to be improving the general reaction 40 There remain several parameters which could be added to the simulation to better mimic the milled system. In addition to the nano-scale mixing (vi), the role of defects created by plastic deformation (ii) could be investigated. Point defects and grain boundaries might enhance the general diffusivity of the system, whereas dislocations could act as nucleation sites 41 . Nanoscale porosity could also change the process dynamics and nucleation features 17 . The combined consequences of all these aspects milling need, however, to be thoroughly studied.

  factors. Crone et al. and Cherukara et al. have indirectly studied activating factors (iv)

L1 2 -FIG. 1 .

 21 FIG. 1. Geometry of the system. (a) Reference system (b) Mixing layer. Dark green = f cc-Ni, light green = f cc-Al, pink = unk-Ni, white =unk-Al where unk stands for unknown atoms. (For the interpretation of the references to color in this figure, the reader is referred to the website where this article appears.)

FIG. 2 .

 2 FIG. 2. Pair correlation function g(r)of the initial configurations after the procedure, normalized on the first peak. g(r) in the case of pure liquid Ni+Al has been added for the sake of information.(For interpretation of the references to color in this figure, the reader is referred to the website where this article appears.)

  . After saturation of the Almelt by Ni-atoms, crystallites of the intermetallic products (NiAl) nucleate at the solidliquid interface and grow along this interface. The microstructure consists of rounded NiAl grains separated by liquid gaps. Further growth of NiAl crystallites occurs due to diffusion through liquid gaps and precipitation at the intermetallic-liquid interface.

Fig. 3

 3 Fig.3tracks the evolution of amorphous-and f cc-atoms in the reference system at 1200 K (REF1200K), above the melting temperature of Al. The rapid increase in the

FIG. 4 .

 4 FIG. 4. (a) Evolution of the Ni concentration in the mixing layer (solid lines), Ni concentration in a 1 nm slice just above or below the interface (dashed lines) and % of B2 formed Ξ in ML20%. (b) Evolution of the number of B2 nuclei (red) and the size of the biggest nucleus (number of atoms) (blue) as a function of time. (c) Snapshot of the mixing layer for REF (top), ML20% (middle) and ML50% (bottom) at 8 ns.

FIG. 5 .

 5 FIG. 5. Snapshots of the atoms at 5 ns for REF and ML20 and at 0.3 ns for ML50. The trajectories of some Ni atoms from the f cc block calculated over the 8 ns simulation are also displayed. Dark green = f cc-Ni, light green = f cc-Al, blue = bcc.

Figure. 5

 5 Figure. 5 shows the diffusion paths during the 8 ns simulation of the Ni coming from the f cc-block. In the case of the REF and the ML20% systems, the grain boundaries act as channels to bring Ni atoms into the middle of the inner layer. As a result, the 50-50 stoichiometry is reached and B2-NiAl continues to grow towards the middle. In the case of ML50%, there is almost no diffusion from the Ni blocks to the middle of the layer, confirming the results presented in Fig. 4a. We next investigated the influence of Ni content in the mixing layer on the mechanism of B2-NiAl formation.

Figure. 6a presents

  Figure. 6a presents the evolution of B2-NiAl in the different systems from 10% to 60%. The inflexion point corresponds to a characteristic time associated with ignition.

FIG. 6 .

 6 FIG. 6. (a) Evolution of at.% of formed B2-NiAl at various Ni concentrations in the mixing layer. (b) Evolution of the initial potential energy of the Al atoms at various Ni concentrations in the mixing layer. The three black vertical lines correspond to the potential energies of the reference structures.

FIG. 7 .

 7 FIG. 7. (a) Evolution of the % of f cc-Al atoms (solid lines) and amorphous Al atoms (dashed lines) in the REF system at different temperatures. Evolution of the at.% of formed B2-NiAl as a function of time at different temperatures. : (b) REF system (c) ML50%.

FIG. 8 .

 8 FIG. 8. Pictures of ML20 at 1200K (top) and ML50 at 900K (bottom).

Figure. 8

 8 shows representative snapshots of the ML20% and ML50% systems at different times during B2-NiAl formation. In ML20%, the seeds nucleate at the interface. The nuclei grow due to precipitation of the Ni and Al atoms present in the inner layer. In the very beginning, some grains may coalesce while others grow, retaining their initial orientation. For an interface oriented normal to the [001] direc-

Fig. 9

 9 FIG. 9. Snapshots of the atoms in the f cc configuration in the inner layer for ML50% at 900K and 1200K and ML60% 1200K. Snapshots correspond to the sate of the system at 8 ns. (Top) Dark green atoms are f cc-Ni atoms and light green atoms are f cc-Al atoms calculated with adaptative-CNA. (Bottom) Local chemical order calculated with Polyhedral Template Matching.

FIG. 10 .

 10 FIG. 10. (a) Evolution of the at.% of B2-NiAl at low temperature (900 K). (b) Evolution of the number of B2 nuclei (red) and size of the biggest nucleus (number of atoms) (blue) over time.

Figure. 10a shows

  Figure. 10a shows the evolution of Ξ as a function of time in the three systems. We noted an unexpected finding: the formation of B2 in the case of the ML40% system is faster and more plentiful than in the two other cases. Figure.10b, in which we note a high nucleation rate, confirms this observation. Grain size then increases due to coalescence and precipitation. This behavior is quite similar to that observed in the 1200K cases.

  Cherukara et al., Manukyan et al. and Xiong et al. studied the shock between pow-

  Kuhn et al. and Ovrutsky et al. have studied nucleation in an Ni 50 Al 50 undercooled amorphous zone

TABLE I .

 I Physical properties of pure elements and intermetallics calculated using Purja Pun and Mishin31 . For more details see Turlo et al.32 .

	Phase	Thermodyamic properties		Structural properties
		T m (K) T m exp. (K)	Cohesive energy (eV) Lattice parameter (nm)
	Ni	1710	1728	-4.45	0.352
	Al	1055	933	-3.36	0.405
	B2-NiAl	1675	1850	-4.51	0.283

TABLE II .

 II Set of simulated systems.

	Generic name	Ni at% a	Temperature (K)	Name
	REF (reference system)	0%	600 K	REF600K
			800 K	REF800K
			900 K	REF900K
			1200 K	REF1200K
	ML (mixing layer)	20%	900 K	ML20%900K
		20%	1200 K	ML20%1200K
	ML (mixing layer)	30%	900 K	ML30%900K
		30%	1200 K	ML30%1200K
	ML (mixing layer)	40%	900 K	ML40%900K
		40%	1200 K	ML40%1200K
	ML (mixing layer)	50%	600 K	ML50%600K
		50%	800 K	ML50%800K
		50%	900 K	ML50%900K
		50%	1200 K	ML50%1200K
	ML (mixing layer)	60%	1200 K	ML60%1200K

a in the mixing layer.

TABLE III .

 III Data on the f cc solid solution formed in different simulations.

	System	ML50% 900K ML50% 1200K ML60% 1200K L1 2 -Ni 3 Al
	Al potential energy (eV)	-4.43	-4.57	-4.62	-5.16
	Ni at. %	50%	55%	62%	75%

  11. The nuclei are labeled 'HON' if they form by homogeneous nucleation and 'HEN' if formed by heterogeneous nucleation. The data presented in Fig.11result from the tracking of several nuclei. The mean square displacement of selected atoms as a function of time is reported in Fig.11a. No difference was noted between Ni and Al atoms. Two types of behavior were observed: (1) A strong and linear increase until 1.5ns followed by a plateau. This behavior is shared by HON40 and HEN30 grains. It reflects high mobility in the system. (2) A slight increase corresponding to very low mobility. This kind of evolution is shared by HEN40 and HON50. The ML40% system was thus found to present both mechanisms of nucleation. Although mobility is hindered close to the interface, this is precisely where heterogeneous nucleation takes place. When mobility is high far from the interface, homogeneous nucleation develops spontaneously. Figure.

ACKNOWLEDGMENTS

The use of computational facilities at the Computing Center of the University of Bourgogne, PSIUN-CCUB, is gratefully acknowledged.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.