A Hémeryck

M Petrantoni

A Estève

C Rossi

M Djafari Rouhani

G Landa

D Estève

Anne Hemeryck
email: hemeryck@laas.fr

A mesoscopic model of the intermixing during nanoenergetic materials processing

Keywords: Multilayers, thin films, interfaces, nanostructures, intermetallic compounds, diffusion

A general mesoscopic model for the simulation of thin film vapor deposition applied to energetic materials, specifically bimetallic multilayers, is presented. We describe the set up of this mesoscopic simulator developed in the frame of a multiscale study by implementing ab initio data into a set of differential equations. We present numerical results relative to the formation of barrier layers as a result of interdiffusion between successive bimetallic AlNi multilayers. The key role of the vacancies species created during deposition is highlighted.

distinct materials [1][2][3][4][5][6][7][8][9] such as Al/Ni, which release heat during mixing. With the thin film vapor deposition techniques, such as Physical Vapor Deposition (PVD) [4,6,10] or Atomic Layer Deposition (ALD) [11], intimate and nearly perfect layers can be processed leading to highly reactive materials [12].

The performances of these energetic materials in terms of stability and reactivity are dependent on the composition, atomic structure, thickness and intimacy of the nanoscale layers. During deposition processes, atomic migrations give rise to an interfacial premixing layer between the two pure polycrystalline constituents [5]. This intermixed layer, in which both composition and structure are unknown, is shown to be one of the major sources of the overall loss of performance of the energetic materials. For instance, it is shown that the intermixing region has a major impact on the thermodynamic driving potential such as the degradation of the self propagating reactions and the reduction of the stored chemical energy [5,10,13,14]. Controlling this interface structure and its formation become crucial in order to get tailored energetic materials in terms of ignition rate and amount of released energy. Our aim in this work is to investigate how the formation of the intermixing layer can be shaped by adjusting the parameters of the deposition process.

In the following, we describe a novel mesoscopic modeling frame in which multilayer deposition is achieved allowing explicit treatment of exothermic reaction with associated dissipation and combustion phenomena. In more details, this mesoscopic model has three main characteristics which make it apart from the other existing models. First, the interfacial layer is not "hand-built" but is formed during the simulation of the deposition process. This simulator could thus be used as a pre-calculation step for a global multi-scale simulation procedure. The objective is to guide other existing approaches [15,16] of the combustion

A c c e p t e d m a n u s c r i p t

4 process based on artificial "hand-built" premixing layer, in order to get a more representative interfacial premixing zone. Second, this model is set up by considering a bottom-up approach [17,18] where ab initio calculations are performed in order to characterize the elementary atomic scale mechanisms (structurally and energetically) occurring during the deposition process. These mechanisms are then used as input parameters in our mesoscopic model. In this manner, the model does not rely on adjustable nor empirical parameters and is only based on theoretical considerations. For sake of brevity, the ab initio data will not be discussed in the present paper. Third, the model is original because it includes mechanisms which are environment dependent with a meso/macroscale model: each atomic reaction is driven by the local composition of the layer in which it occurs. This predictive model is able to establish a direct link between the species concentrations in each layer of the material and their relation with the composition of other layers and the processing parameters. Although the model is generic for bilayer materials, we detail, in this paper, a set of equations intended to be validated on the model system composed of a bimetallic material, i.e. an AlNi multilayer film.

Methodology: Mesoscopic modeling

Atomic scale mechanisms

Using ab initio calculations, we have identified eight general elementary reactions that can occur in the multilayered Al/Ni materials: three of them are described as intra-layer reactions and the five others are inter-layer reactions. On the other hand, these reactions are qualified as: (i) "creation" like a defect formation mechanism, (ii) "exchange" between two chemically distinct atoms and (iii) "diffusion" mechanisms through the layers. We denote respectively as SAl and SNi, aluminum and nickel atoms positioned on a lattice site of the fcc crystals; as IAl and INi, aluminum and nickel atoms in a non-lattice position, and as V a vacancy. This vacancy can be initially present in the layer, created by atomic migration between the layers or provided to the system during the deposition process.

A c c e p t e d m a n u s c r i p t 5

We detail in the following the intra-layer reactions occurring in a layer referred to as n , and inter-layer reactions occurring between two adjacent layers referred to as n and 1 r n . These elementary mechanisms are written as:

-Vacancy formation [19,20], providing a vacancy (V) and a non-lattice atoms (IAl or INi) [21] from an atom in a lattice position (SAl or SNi), is considered for both constituents:

V(n) INi(n) SNi(n) l (R.I) V(n) IAl(n) SAl(n) l (R.II)
This reaction is considered as possible in the model despite its large formation energy for Ni-Al system, estimated at roughly 8 eV using ab initio calculations, since the model is generic and does not refer to this special case. The occurrence of this reaction is however automatically inhibited in Ni-Al system because of its very low reaction rate.

-Exchange between two kinds of chemically different species, respectively in lattice and non-lattice positions

SAl(n) INi(n) IAl(n) SNi(n) l (R.III)
-Exchange of two crystalline species present in two adjacent layers is allowed:

SAl(n) 1) SNi(n 1) SAl(n SNi(n) r l r (R.IV)
-Two diffusion mechanisms mediated by the presence of a vacancy

V(n) 1) SNi(n 1) V(n SNi(n) r l r (R.V) V(n) 1) SAl(n 1) V(n SAl(n) r l r (R.VI)
-Two mechanisms relative to the diffusion of non-lattice species

1) INi(n INi(n) r l (R.VII) 1) IAl(n IAl(n) r l (R.VIII)

Implementation of the Chemical Rate theory

In order to access to the temporal evolution of the system during the technological process, it is necessary to follow the concentration of each species constituting the layer. The species concentrations are noted as > @ SNi , > @ SAl , > @ INi , > @ IAl and > @

V for the five identified species.

Two kinds of mechanisms have been distinguished:

-The "incoming molecular flux" providing species during deposition and filling in the -For all others implemented mechanisms, the performed ab initio calculations reveal that an activation energy is required. This energetic barrier determines the probability of occurrence of a given event that depends also on the temperature. This is implemented in the model using the Arrhenius law which leads to the reaction constants i k :

layer g i T R A i e ȣ k u With i
A the activation energy of the mechanism i (J.mol -1 ), g R the ideal gas constant (8.31 J.mol -1 .K -1 ) and layer T the time dependent layer temperature (K). ȣ is a dynamic prefactor corresponding to the Debye frequency. ȣ represents an attempt frequency of jump and is estimated at 10 13 s -1 .

The effects of phase changes in the film such as alloy formation due to the melting of the reactants or products formation are neglected.

Implementation of heat exchange

The chemical reactions occurring in the energetic materials are highly exothermic and produce heat during mixing. This heat release is also present during the deposition process due to atomic migration leading to a structural rearrangement and to the formation of the premixing zone. In order to build a realistic and accurate model dedicated to the energetic materials, the temperature variations in the multilayer films have to be taken into account.

We assume a uniform time dependent temperature throughout the multilayer film. This approximation is mainly justified by the small thickness of the multilayer films, of the order of few nm, compared to the thickness of the silicon substrate, of the order of 100 Pm. Spatial variations of the temperature, if present, occur only in the silicon substrate. They could be taken into account via heat conduction in the silicon substrate film. However, these spatial Also, in the present version of the model, physical properties such as thermal conductivity, density and heat capacity are assumed to be independent of the temperature and of the layer composition. In a future version of the model, the temperature dependence of these parameters can be easily taken into account, if necessary. Considering the above points, the temporal evolution of the temperature can be written as:

) T ı7 ) T h(T Q dt dT e c ȡ 4 0 4 0 i i Si Si Si ¦
With Si ȡ the silicon density (2330 kg.m -3 ), Si c the silicon heat capacity (700 J.kg -1 .K -1 ), h the convective heat exchange coefficient, Si e the silicon substrate thickness (m), Si T is the temperature of the silicon substrate, ı the Stefan-Boltzmann's constant (5.67x10 -8 W.m -2 .K - 4 ), 0 T the ambient temperature and i Q the heat released by each exothermic reaction during the time interval dt. This last parameter is environment dependent (Section 2.4).

The model is therefore composed of 1 5n differential equations for the description of the concentration variations for a material with n layers: five equations for the five species in each layer and one equation for the uniform temperature.

Implementation of the environment dependent reactions

The originality of the model is that the activation barriers i A are environment dependent: they explicitly depend on the composition of each layer and account for the heat production

A c c e p t e d m a n u s c r i p t

9 noted as i Q . In order to implement this environmental dependence in the model, the total energy is assumed to be the sum of an electronic term and nearest neighbor pair interactions.

The energy of an atom of type J can then be written as:

jk k jj j 0 İ n İ n (J) E J E
Where 0 E corresponds to the electronic term, j n and k n are the ratios of species j and k in the layer, and jj İ and jk İ are the pair interaction energies between two nearest neighbors, of the same type or of different types, respectively. The ratios j n and k n of species in the layer are expressed as > @ > @ > @ > @ V K J J n j and > @ > @ > @ > @

V K J K n k .
Where [J], [K] and [V] respectively stand for the concentrations of J and K species, and for the vacancy concentration. We have assumed that vacancies may be present in the layers as a result of kinetic effects and the lattice mismatch between the two materials.

The energy values for the starting position, the final configuration and the transition state are evaluated at each step of the simulation, for all possible elementary reactions and for each layer in the film. Using these calculated energies (Figure 1), the activation energies i A as well as the released heat i Q for both the direct and the reverse reactions and for each elementary mechanism can be determined. In this way, the knowledge of only three energetic parameters dependent on the composition of each layer, allows us to address the desired environmental dependence.

Preliminary results

In the following, we present preliminary results concerning the Ni/Al multilayer bimetallic films deposited on silicon substrates. The aim of the simulation is a first validation of the model through a simple system for which experimental and theoretical data are available in the literature. Also, the technology process conditions have been simplified as much as For both cases, deposition of Ni or Al, we consider that the first material is already deposited so that the starting point of the simulation concerns the deposition of the second material. or Ni/Al(111), with a growing multilayer film maintained at ambient temperature, no interfacial mixing can be obtained if no vacancy is provided during the deposition process. In fact, the considered atomic scale reactions are very difficult to occur in these metallic materials at low temperatures, as expected from DFT results.

For instance, if we focus on the exchange between two atoms, SNi or SAl, described in the reaction R.IV, the exchange of SNi located in a pure nickel layer with a SAl located in a pure aluminum layer, will require energy barriers up to 5 eV. This has been determined using ab initio calculations where the theoretical minimum energy pathway is characterized by a multi- In a second set of runs, the simulations are performed by considering a perfect material composed of four Ni monolayers in direct contact with four Al monolayers and a film temperature maintained at 298 K. These simulations run last for 1s with no incoming gas flux.

When an artificial amount of vacancies, 10 % here, is introduced in the interfacial monolayer between the two pure constituents, the mixing reactions are activated, even at room temperature. Here, we wish to emphasize that the vacancy presence in the interfacial layer will be key in the mixing taking place during the vapor deposition as the vacancies favor the mixing at the interface. Indeed, activation barriers, as calculated using the three states (Figure 1), are strongly reduced in the presence of vacancies.

Conclusion

During multilayer deposition processes, an intermixing reaction inevitably occurs between the two constituents leading to the formation of a premixed zone which separates pure layers of the films. This paper presents a first attempt to set up a model for the deposition process and the structural arrangement during the deposition, based on a multi-scale approach combining chemical rate theory and ab initio results. In these preliminary results, a first trend appears:

the diffusion of Ni species seems to be easier across the bulk compared to the Al species diffusion. The key role of vacancies species provided to the system during the deposition process is also highlighted. We reveal that if no vacancy supply is provided during deposition process, for multilayer films maintained at ambient temperature, it is tricky to form an interfacial premixing layer due to the high activation barriers required to allow atomic reactions. We have shown that, in order to develop a realistic model of the deposition process, these developments, the model presented here could be coupled to other existing modeling approaches as they mainly focus on the heat loss phenomena due to the intermixing layer, or on the combustion process. These models, starting usually from empirical premixed layers, can use our approach to define process dependent initial configurations.

List of references

A c c e p t e d m a n u s c r i p t 1 .

 1 IntroductionNanostructured reactive multilayer films are typically formed by alternating chemically

A c c e p t e d m a n u s c r i p t 6 These

 6 reactions have been characterized using ab initio/Density Functional Theory. Five different systems have been considered. Two are pure Al and Ni bulk materials, one is AlNi bulk alloy and the two remaining represent Al(111) and Ni(111) surfaces. Each system consists of a supercell of 96 atoms distributed in 6 layers of 16 atoms. Periodic boundary conditions are applied in all three directions. These results are still the subject of current research and should be published soon.

  empty layers. Ab initio calculations have shown that both atomic adsorptions of Ni on Al(111) and Al on Ni(111) result in hexagonal lattice positions SNi and SAl. The adsorbed atom is placed above adjacent sub-layer bulk lattice atom. The incoming molecular flux is distributed among successive deposited layers, proportionally to the area exhibited by each layer. The adsorption mechanism which does not require any activation energy is driven by the Maxwell Boltzmann statistics in the gas phase, pressure of the considered gas (Pa), S the surface of the elementary cell (considering Al or Ni lattice parameter) (m²), gasM the molar mass of the gas (kg.mol -1 ), gas T the gas temperature (K), and C a constant. In this manner, the pressure and the temperature of the incoming gas are introduced as in a real processing experiment.

A c c e p t e d m a n u s c r i p t 7

 7

A c c e p t e d m a n u s c r i p t 8 variations

 8 are not implemented in this first version of the model, since we consider only mechanisms where the atomic species diffusivities are smaller than diffusivity. In this model, temperature variations result from exothermic reactions and the convective and radiative heat losses on the wafer surfaces. Concerning the radiative heat losses, the multilayer film is approximated as a black-body radiating with unit emissivity, which ideally exchanges with unbounded surroundings at ambient temperature.

A c c e p t e d m a n u s c r i p t 10 possible.

 10 For instance, the temperature has not only been considered as uniform over the whole sample, but also constant in time and equal to the ambient temperature. This approximation is justified in cases where the silicon substrate is strongly cooled during metal deposition.In all following simulations, the ambient temperature is kept fixed at 298 K. The incoming gas is at 350 K under a pressure of 0.005 Pa. We mainly focus on the role of vacancies present at the interface between the two nanolayers. Our objective is to understand to what extent their presence is necessary for the formation of an intermixed zone between the pure nanolayers. In a first set of simulation runs, we successively investigate the deposition of Ni/perfect Al(111) substrates and of Al/ perfect Ni(111) substrates, to get insight on both interfaces. However, in these preliminary simulations, we have not examined multilayer deposition. All simulations of this set run for 3s with μm 1 e Si .The atomic pair interaction energies have been estimated via ab initio calculations. Five periodic structures of pure Al, of pure Ni and of AlNi alloy have been considered using a large spectrum of situations: vacancies, interstitials, Frenkel pairs, atom/surface interactions... The values for the electronic energies and the pair interaction energies have been estimated by fitting the total energies to the corresponding values obtained using DFT, in all considered cases. The best fit is observed for: lattice positions (INi-INi, IAl-IAl, IAl-INi) have not been considered in this first version of the model, due to the low probability of finding two interstitials close to each other.

FigureA c c e p t e d m a n u s c r i p t 11 layers

 11 Figure 2a represents a schematic view of the starting configuration and the partially covered

A c c e p t e d m a n u s c r i p t 12 Al

 12 step reaction involving several non-lattice positions for the migrating species and its neighboring atoms. Despite the small intermixing occurring during room temperature simulations, we can still observe two general trends, in agreement with ab initio calculations performed to determine the adsorption pathways. The first concerns the heat loss during deposition which is larger when Al is deposited on a Ni substrate. The second is that the intermixing is asymmetric: Ni atoms penetrate easily the Al(111) substrate, while Al remains on Ni(111) surface. The theoretical pathways obtained show that only a 0.06 eV activation barrier allows the penetration of the Ni atom into the first Al(111) sub-layer. The activation barrier for the migration of the Ni atom into the second Al(111) sub-layer is of 1.2 eV, in the range of bulk migration energies. Looking at the Al on Ni(111) case, the activation barrier for the Al atom to penetrate the first Ni(111) sub-layer amounts to 1.56 eV, much higher than the previous value. This reaction is practically forbidden at room temperature. Therefore, the intermixed zone, for room temperature processes, results at the Ni/Al(111) interface rather than at the /Ni(111) interface. This asymmetry shows that tuning experimental conditions is a key for the optimization of the multilayer films, particularly in terms of formation of intermixed zones. These conditions might eventually be different for the two, Ni and Al, deposition processes.

A c c e p t e d m a n u s c r i p t 13 it

 13 is necessary to investigate and implement the mechanisms of supply of vacancies to the interfacial region. As a result of the preliminary results discussed in this paper, future work includes an accurate and realistic rate of vacancies supply during deposition process. Beyond

A c c e p t e d m a n u s c r i p t 15 Figures caption Figure 1 :

 15caption1 Figures caption

Figure 2 :A c c e p t e d m a n u s c r i p t 16 Figure 1 A c c e p t e d m a n u s c r i p t 17 Figure 2 A c c e p t e d m a n u s c r i p tFigure

 2161172 Figure 2: (a) Schematic views of the starting configurations of the multilayer films before

Table of the length of the manuscript

	Words contained in the text (abstract + list of	3539 words
	ref but not figure caption)
	Number of figure


A Perron

S Garruchet

O Politano
email: politano@u-bourgogne.fr

G Aral

V Vignal

Oxidation of nanocrystalline Aluminum by variable charge molecular dynamics

Keywords: A, metals A, oxides A, thin films D, diffusion D, microstructure

We investigate the oxidation of nanocrystalline aluminum surfaces using molecular dynamics (MD) simulations with the variable charge model that allows charge dynamically transfer among atoms. The interaction potential between atoms is described by the electrostatic plus (Es+) potential model, which is composed of an embedded atom method potential and an electrostatic term. The simulations were performed from 300 to 750K on polycrystalline samples with a mean grain size of 5 nanometers. We mainly focused on the effect of the temperature parameter on the oxidation kinetic. The results show that, beyond a first linear regime, the kinetics follow a direct logarithmic law (governed by diffusion process) and tend to a limiting value corresponding to a thickness of ~3nm. We also characterized at 600K the effects of an external applied strain on the microstructure and the chemical composition of the oxide films formed at the surface. In particular, we obtained a partially crystalline oxide films for all temperature and we noticed a strong correlation between the degree of crystallinity of the oxide film and the oxidation temperature.

Introduction

Metal/oxide interfaces play critical roles in many applications, such as structural or functional materials in materials science, microelectronics, chemical and biomedical applications. In materials science for instance, the following technological materials are based on metal/oxide interfaces: functional ceramics with metals (cermets), oxide dispersion-strengthened alloys, oxide coatings on metals functioning as thermal barriers or natural corrosion protection layers. This is in particular the case of oxide passive films formed spontaneously on aluminum and aluminum alloys which protect the surface against further oxidation and corrosion. Also, because of their specific physical properties (i.e. large dielectric constant and height barrier for electron tunneling), aluminum oxides are promising candidate for gate dielectrics in the nanoscale CMOS technology, magnetic tunneling junctions, magnetic sensors, magnetic storage devices and solid-state memory, etc. [1][2][3]. At the same time, significant progress has been made to understand the combustion of nano-aluminum powders and the highly exothermic reaction of fine-grained aluminum particles (size < 100 nm) with oxygen because of their use as potential additives to propellants for rockets [4][5][6][7][8][9]. For all these applications, it is crucial to understand the physico-chemical processes which control the reactivity of the metallic substrate and the thickness, morphology, chemical composition and microstructure of the oxide film [1].

The development of numerical studies of metal/oxides interfaces as well as oxide growth on the top of metals have required the development of the new class of interaction potential such as for example the variable charge method proposed by Streitz and Mintmire which is often referred as the (Es+) potential [10]. This last combines an embedded atom method (EAM) potential and an electrostatic term to model the metal-metal, oxygen-oxygen and metaloxygen interactions. The EAM contribution is predominant in the metallic region whereas the which is the thermodynamically stable Al oxide. Later, Streitz and Mintmire also showed that the (ES+) potential was able to describe the gamma alumina (γ-Al 2 O 3 ) structure with a reasonable degree of accuracy [11]. This last point is of importance in the framework of the present study as high temperature oxidation of aluminum usually leads to the formation of transition alumina depending on the temperature.

Initially designed to simulate the properties of one metal (namely Al) interacting with oxygen, the Streitz and Mintmire potential was recently extended/modified to simulate many metals (Ti, Zr, Fe, Ni, Co,…) and alloys [12][13][14][15][16][17][18]. It was used to study nanocrystalline metal-oxide properties [16][17][18] or processes involved in the oxidation of metallic substrates [12,13,[19][20][21][22].

In this last case, the capability of the interatomic potential to determine the atomic charges on atoms due to its local environment was used to model the chemical reactions (bond formation and bond breakage) between a metal and an oxygen gas within a molecular dynamics simulation.

Recently, a new approach [23] was developed to solve one deficiency of the afore-mentioned variable charge models which were found to weakly describe the iono-covalent nature of the metal-oxygen bonding. For example, they lead to the underestimation by 30% of the calculated formation energy of an oxygen vacancy. This new generation of variable charge models incorporates new physics through a tight binding analytical model [24] and gave very satisfactory results for ZrO 2 and TiO 2 polymorphs. It looks very promising for the development of the future simulation of metal/oxyde systems.

A c c e p t e d m a n u s c r i p t

Oxidation of Al nanoclusters [19,20] or single-crystal surfaces of Al [21,22] for pressure ranging from 10 to 40 times the normal state were studied by variable charge molecular dynamics (VCMD) approaches. Under their oxidation conditions, these works clearly characterized the formation of an amorphous oxide layer. They reported a growth mechanism proceeding by both anionic and cationic diffusion and a direct logarithm relationship for the oxidation kinetic. Note that these VCMD simulations never permit to characterize the formation of a well crystalline oxide film.

On the other hand, experimental investigations on the high temperature oxidation of pure In the present approach, we studied the oxidation under a constant molecular oxygen density (1 atm) of a polycrystalline Al substrate with a mean grain size of 5nm. We focused on the relation between the growth kinetics and the evolution of the chemical composition, by studying the structure of the developing oxide scale. In order to compare with the existing

A c c e p t e d m a n u s c r i p t

experimental results of Jeurgens et al., we used the same range of oxidation temperature (i.e., 300, 600 and 750K). Moreover, to highlight the influence of an applied strain on the oxidation mechanisms, one atomistic system was oxidized at 600K after an uniaxial deformation (+7%).

Computational Methods

We used the variable charge molecular dynamics method (VCMD) originally proposed by

Streitz and Mintmire for Al-O systems [10]. In this formalism, the interatomic potential is a sum of the electrostatic and Finnis-Sinclair embedded-atom potentials and allows the atomic charges to vary dynamically in response to changes in the local environment of the atoms.

The variable-charge approach determines atomic charges at each molecular dynamic step in such a way that the electrostatic potential energy is minimized, with the charge-neutrality constraint, i.e. the sum of the atomic charges is equal to zero. This condition is equivalent to the electronegativity equalization principle which requires that each atoms, aluminum (Al)

and oxygen (O), have an equal chemical potential.

The VCMD code used in this approach is developed by the "Collaboratory for Advanced

Computing and Simulations" group at the University of Southern California [20]. It is highly parallel and allows to simulate million atoms metal-oxide systems. To speedup calculations, state of the art techniques are implemented: multiple time steps, multilevel preconditioned conjugate gradient method, fast multipole method and domain decomposition approach are used to implement the MD algorithm on parallel architectures [30][31][32].

The initial Al polycrystalline sample is a cube of size 10.256 nm containing 62 932 atoms.

The bulk polycrystalline system was created at 0K by using a procedure analogous to the one

A c c e p t e d m a n u s c r i p t

used to obtain Voronoi cells (Fig. 1a). It was then thermalized by successive jumps of 50K until the desired temperature by using Noose-Hoover thermostat chain [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. The slab geometry is then created by incorporating in the simulation box two vacuum slabs along the Z-direction while periodic boundary conditions are imposed along X and Y axis. For more details about the sample preparation procedure, the reader should refer to [START_REF] Perron | [END_REF]. The final dimension of the simulation box is Z = 15nm and X = Y = 10.256nm (Fig. 1b). Our MD simulations were performed by using a constant oxygen molecular density equal to the normal state P 0 (1atm) above the Al polycrystalline surface. Every MD steps, the number of oxygen atoms present in the gas is examined and if this number is less than the imposed value, then new oxygen molecules (with a velocity equal to the oxidation temperature) are randomly introduced into the region above the Al surface in order to maintain the molecular oxygen density constant.

During the oxidation stage, only the Al atoms belonging to the lower 0.7nm part of the polycrystalline substrate (Fig. 1b) were kept thermalized by using the thermostat. The constrain on the other Al atoms was removed and their kinetic energy (temperature) was free to evolve during all the oxidation process [START_REF] Perron | Apport de la dynamique moléculaire et de techniques expérimentales à une échelle locale pour l'étude de surfaces d'aluminium déformées et/ou oxydées[END_REF].

We integrated the equation of motions by using a multiple time step algorithm in which the MD time step is dt = 2fs for short-range forces and dt = 5fs for long-range forces. As the calculation of the charges is the most consuming part of VCMD, the local atomic charges are only calculated every 100th time steps. That is, after their updates the charges are kept constants during the following 99 steps to compute the forces and energies. A cutoff radius of 0.6nm was used to compute the Finnis-Sinclair potential.

A c c e p t e d m a n u s c r i p t

The strain samples were created by following a procedure which gave satisfactory results to study the multilayer surface relaxation of Al single-crystals [START_REF] Politano | [END_REF] or the formation of surface roughness on strained polycrystalline Al substrates [START_REF] Perron | [END_REF]. It consists of a direct deformation of the sample to the desire length. That is, we rescaled instantaneously the positions of the atoms along the Y-direction in order to obtain the required strain. The sample was then relaxed during 50 ps.

Results and discussions

Kinetic law

The kinetic curves obtained for all systems (unstrained : 300, 600 and 750K ; strained : 600K) are plotted on Fig. 2. In this figure, the y-axis represents the evolution of the number of oxygen atoms introduced into the simulation box per unit surface area (nm 2 ). It can be assimilated to the variation of sample weight measured experimentally (Fig. 2a). The thickness of the oxide film is determined by calculating the difference between the mean position in the Z-direction of the ten highest Al atoms (near the gas/oxide interface) and the ten deepest oxygen atoms (at the metal/oxide interface). As one can notice, Fig. 2a and Fig. 2b are similar to the classical representations obtained during high temperature oxidation experiments: sample weight gain/loss and film thickness evolution.

The kinetic curves presented on Fig. 2a show and identical behaviors for all oxidizing conditions. They can be divided into two parts: a first linear law followed by a direct logarithmic relationship. According to the oxidation theory and heterogeneous kinetics [START_REF] Birks | Introduction to the high-temperature oxidation of metals[END_REF][START_REF] Sarrazin | Mechanisms of High Temperature Corrosion: A Kinetic Approach[END_REF][START_REF] Khanna | Introduction to High Temperature Oxidation and Corrosion[END_REF][START_REF] Soustelle | Cinétique hétérogène : Mécanismes et réactivité[END_REF], during the linear part the interfacial processes (dissociation of oxygen molecules, sorption of oxygen atoms, Al-O reaction) are predominant. Indeed, at the beginning of the oxidation, the oxygen atoms react directly with the Al substrate and the rate of this reaction

A c c e p t e d m a n u s c r i p t

governs the global kinetic. As presented on Fig. 2a, the linear law is available until around 1ns of oxidation and corresponds to an oxide thickness of ~1.5nm (Fig. 2b). After this point, the oxygen atoms do not have direct access to the Al substrate but have first to cross the oxide film which acts as a barrier. So, this second kinetic part of the oxidation is controlled by diffusion of atomic species (oxygen and aluminum) through the film.

As two different processes occur during oxidation, we analyzed the kinetic curves by using a comprehensive law combining linear and logarithmic terms:

N/S = K l t + K log log(D t + E) .
Eq. 1

In this expression, t is the oxidation time, N/S is the number of oxygen atoms per unit area in the oxide film, K l and K log respectively the linear and logarithmic constants. D and E are mathematical parameters. Such combination of kinetic laws is rather standard when different processes occur during oxidation (the best example is the well known linear-parabolic relationship) and is detailed in textbooks on high temperature oxidations [START_REF] Birks | Introduction to the high-temperature oxidation of metals[END_REF][START_REF] Sarrazin | Mechanisms of High Temperature Corrosion: A Kinetic Approach[END_REF][START_REF] Khanna | Introduction to High Temperature Oxidation and Corrosion[END_REF][START_REF] Soustelle | Cinétique hétérogène : Mécanismes et réactivité[END_REF].

The numerical fitting of the kinetic curves by using Eq. 1 is reported on Fig. 3. We can observe an excellent fit of the oxidation kinetic by the comprehensive law: for low values of N, the linear law is predominant and then, with increasing N, the logarithmic law becomes majority. The parameters obtained from the numerical fitting are given in Table I. As one can notice, K l which is the slope of the linear part remains almost constant for all oxidizing conditions. Therefore, sorption phenomena of oxygen are not modified by oxidation temperature or applied strain. On the other hand, K log increases with temperature and applied strain. Therefore, the diffusion of atomic species is facilitated by increasing the motion of atoms (temperature) and the inter-atomic distances (deformation). The activation energy, which can be deduced from the Arrhenius plots (ln(K log ) versus 1/T), was estimated at 1.3kJ.mol -1 . This value is rather low compared to experimental results ranging from 9.5 to

A c c e p t e d m a n u s c r i p t

238kJ.mol -1 for the oxidation of bulk Al or nano-particles (~50nm) [7][8][9]. This difference might be explained by the very high reactivity of our samples induced by their small nanocrystalline microstructure (5nm grains), their total purity and the absence of a protective alumina scale at the surface of the sample. Finally, we note that the growth processes proceed almost exclusively by anionic migration (internal movement of oxygen anions through the oxide to the metal substrate).

The presence of the two kinetic regimes is corroborated by previous studies from Hasnaoui et al. on the oxidation of single-crystal Al surfaces. Indeed, these authors mentioned a short "transition" initial stage before the direct logarithmic law [21,22]. Based on the present work,

we are now able to conclude that the "transition" stage observed by Hasnaoui et al. was linear and governed by interfacial processes.

Concerning the oxide film thickness represented on Fig. 2b, we can observe an increase with the temperature until ~3nm. This agrees well with the value of 5nm reported by Jeurgens et al. with in-situ XPS and Auger spectroscopy [27]. Furthermore, small variations of thickness are observed after 1ns for 600 and 750K. We observed a layer by layer growth mode for the oxide film at 300K. At higher temperature, the growth of the oxide is less uniform and leads to the fluctuations observed on Fig. 2b. Furthermore, microstructural changes of the oxide film can also produce this behavior. This last will be developed in the next section.

Chemical composition

As suggested above, the morphology and composition of the oxide films seems to vary with the temperature. To clarify the changes in composition, the fractions of aluminum, oxygen atoms and the various atomic densities have been calculated at the end of each oxidation. In order to avoid any spurious effect due to the interfaces which are too

A c c e p t e d m a n u s c r i p t

concentrated in aluminum or oxygen atoms, the calculations are performed between from 0.5nm under the oxide/gas limit to 0.5nm above the substrate/oxide limit. It corresponds to the middle part of the oxide film. The obtained values are reported in Table II and are discussed below.

Concerning the temperature: an increase of the aluminum atomic fraction in the oxide film is observed between 300 and 600K (and consequently a proportional decrease of the oxygen fraction). Moreover, the Al atomic density decreases slightly with temperature whereas the O atomic density decreases strongly. As a consequence, the total density is decreasing between 300 and 600K. To resume, the increase of the oxidation temperature induces an increase of the kinetic law and produces oxide films which are less dense and enriched in Al. Concerning the oxide film formed at 750K, the composition is similar to the one formed at 600K but the density is different. It indicates that the microstructure changed between 600 and 750K. This point will be studied in more detail afterwards.

Concerning the applied strain (+7%): we noticed a slight increase of the O fraction and a strong increase of the O atomic density in the oxide film. To our point of view, these observations are due to the increase of the inter-atomic distances (Al-Al) which facilitates the insertion of O atoms in the substrate.

Microstructure characterizations

The temperature of the oxide film was computed as a function of time and is shown on Fig. 4. First, we observed a slight decrease of the temperature until the oxide scale cover completely the surface. After this very short transition regime, a great increase of temperature in the oxide is observed until ~1.5ns. This phenomenon corresponds to the linear kinetic law determined previously (sorption process). After 1.5ns, which corresponds to the beginning of

A c c e p t e d m a n u s c r i p t

the logarithmic kinetic law (diffusion processes), the temperature of the oxide films reaches a maximum and then decreases very slowly with time. So, the evolution of the oxide temperature can be associated to the number of oxygen consumed. The oxidation reaction of Al is highly exothermic (-1 675kJ.mol -1 ) and, therefore, leads to a great increase of the oxide film temperature (included between 1 000 and 2 000K). These high temperatures will favor the formation of crystalline Al 2 O 3 structures as it increases the atomic motion/rearrangement in the films.

Pair distribution functions (PDF) for the Al-O bonds have been calculated in the middle part of the oxide films in order to characterize their microstructure. As presented in Fig. 5, many peaks are observed on the PDF and indicate the presence of long range order (a partially crystalline structure). But, no clear difference can be seen between the different curves.

Moreover, one cannot determine the crystalline microstructure of the films from Fig. 5. A second step consists in calculating the environment of each atom by applying a common neighbor analysis on Al and O atoms [START_REF] Ackland | [END_REF]42]. The crystallographic structure of the thermodynamically stable alumina (α-Al 2 O 3 ) is based on an hexagonal network of O and Al atoms [43]. Therefore, as a first approach we used the percentage of oxygen atoms located in hexagonal compact position (hcp) within the oxide as a criteria to characterize the cristallinity of the film formed at the Al surface. The results are reported Table III. As one can see, the percentage of O in hcp positions increases with temperature and also by applying a strain on the metallic substrate. This last highlights an increase of the degree of crystallinity from ~32% at 300K to ~43% at 750K. Finally, the change in microstructure between 600 and 750K suggested previously is confirmed. However, many transition alumina exist -γ, η (cubic), θ, θ', θ'', λ (monoclinic), δ-Al 2 O 3 (tetragonal or orthorhombic), α (trigonal), κ (orthorhombic)

, χ (hexagonal) -and the present study did not allowed us to determine which alumina polymorphs is present. In a forthcoming study, a more extended quantification of the ratios of Concerning the polycrystalline microstructure, we did not observed a significant difference of reactivity between the Al atoms located in grains or in grains boundaries. Also, no preferential diffusion of the oxygen ions was observed at the grain boundaries.

Conclusion

We performed VCMD simulations on the oxidation of nanocrystalline aluminum between 300 and 750K. The influence of an external applied strain was also studied at 600K. This last permit us to obtain the corresponding kinetic curves and the formation of an thin oxide film which was analyzed with a common neighbor analysis and with pair distribution functions.

The kinetic curves show the formation of a film with a limited thickness of about 3nm. Those curves can be divided in two regimes: a linear and a logarithmic. According to heterogeneous kinetics, these two regimes are governed by interfacial processes and by diffusion respectively. A comprehensive kinetic relationship permit us to compute the linear and a logarithmic constants. The linear constant was found to be independent of the temperature and the applied strain. At the opposite, the logarithmic constant, related motion of the atoms in the oxide scale, increases with both temperature and strain. Also, an increased of the degree of crystallinity within the oxide was observed with temperature and strains. We did not observed any preferential diffusion at the grains boundaries in the polycrystalline substrate.

Moreover, it is important to note that previous MD studies on single-crystal aluminum samples only reported the formation of amorphous alumina films. Here, we obtained the formation of a partially crystalline structure (42% of O atoms in hcp positions at 750K). To

A c c e p t e d m a n u s c r i p t

our opinion, this might be related the lower oxygen pressure (P = 1atm) used in the present approach as well as to the use of a large surface area which reduce the border effects. Also the kinetic factor is crucial and in order to leave the atoms reorganize to obtain a well crystalline oxide, we will surely have to simulate much longer times.

A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p t

  coulombic potential is used to model the ionic interactions within the oxide. It was initially fitted to model several properties of metallic aluminum (Al) and alpha alumina (α-Al 2 O 3 )

  aluminum surfaces are difficult to perform due to the very high reactivity of Al with O and the presence of the native protective alumina scale. Surface preparation under high vacuum conditions is required and in-situ experimental techniques have to be conducted under very low oxidation pressure. Many experiments on bulk aluminum surfaces were developed by Jeurgens et al. with XPS and Auger spectroscopy. These studies were performed at very low oxidation pressure (∼10 -4 Pa) and mainly focused on the reaction kinetic, the composition and the structure of the alumina scale as a function of the temperature. That is i) at low temperatures (T < 573 K), an amorphous and Al-deficient (as compared to γ-Al 2 O 3 ) thin oxide film (< 5 nm) is formed; ii) at higher temperatures (T≥673 K), during the initial fast oxidation stage an Al-enriched amorphous oxide film is obtained and during the second slow oxidation stage the films gradually reach the stoechiometric composition and the crystalline structure of γ-Al 2 O 3 [25-29].

A c c e p t e d m a n u s c r i p t tetrahedral AlO 4

 4 and octahedral AlO 6 units will certainly help us to determine which alumina polymorph is formed during the simulation.

Fig. 1 .

 1 Fig. 1. (a) Bulk polycrystalline aluminum sample composed of 62 969 atoms. The grey light

Fig. 2 .

 2 Fig. 2. (a) Evolution of the number of oxygen atoms consumed by the oxidation process and

Fig. 3 .

 3 Fig. 3. Oxidation time as a function of the number of oxygen atoms. The symbols are the

Fig. 4 .

 4 Fig. 4. Temperature of the oxide film formed as a function of time for different oxidation temperature (300, 600 and 750K) and straining conditions (+7%).

Fig. 5 .

 5 Fig. 5. Pair distribution function, g(r), for the Al-O inter-atomic distances in the middle part of the oxide films.

A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p t FIGURE 1 ZFIGURE 3 A c c e p t e d m a n u s c r i p t FIGURE 4 A c c e p t e d m a n u s c r i p t

 134 FIGURE 1

FIGURE 5

 5 FIGURE 5

Table I .

 I Parameters values of the comprehensive law according to the oxidation temperature and applied strain.

	K l	300K	D		K log 600K	E 750K
	atoms.nm -2 .ps -1 0.0634 Applied strain 300K 0% 600K 0.0655 750K 0.0634 Al atomic fraction in the 0.50 oxide film (%) 600K -+7% 0.0592 O atomic fraction in the oxide film (%) 0.50	ps -1 0.0209 0.0197 0.0133 0.0126	0% 0.53 0.47	atoms.nm -2 57.4320 ε X = +7% 62.2611 85.6410 0.51 76.588 0.49	0% 0.00932 0.02868 0.52 0.4500 0.3095 0.48
	A c c e p t e d Aluminum atomic density (atomes.nm -3 ) 49.8 Oxygen atomic density (atomes.nm -3 ) 49.0 Total atomic density (atomes.nm -3 ) 98.8	m a n u s c r i p t 49.1 49.2 43.5 47.3 92.6 96.5	48.5 44.7 93.2

Table II .

 II Evolution of the atomic fractions and atomic densities in the oxide films formed at 300, 600 and 750 K, and after deformation (+7%).

	Conditions	300K	600K	750K	600K -Strained
	Degree of crystallinity	31.7%	37.6%	41.5%	42.5%
	A c c e p t e d	m a n u s c r i p t

Table III .

 III Percentage of oxygen atoms in hexagonal compact position within the oxide layer.

Acknowledgments

The authors want to thanks P. Vashishta, A. Nakano and R. Kalia for fruitful discussions. This work has been partially performed under the HPC-EUROPA (RII3-CT-2003-506079) with the support of the European Community-Research Infrastructure Action under the FP6 "structuring the European Research Area Progamme". Part of computations have been carried out at the CRI-CCUB, IDRIS and ULAKBIM high performance computing center.

A c c e p t e d m a n u s c r i p t