Lu Xie

Pascal Brault
email: pascal.brault@univ-orleans.fr

Jean-Marc Bauchire

Anne-Lise Thomann

Larbi Bedra

Running head: Cluster and thin film growth by MD simulations

Keywords: 81.15.Aa, 81.15.Cd, 68.55.J-, 02.70.Ns

Carrying out molecular Dynamics simulations is a relevant way for understanding growth phenomena at the atomic scale. Initial conditions are defined for reproducing deposition conditions of plasma sputtering experiments. Two case studies are developed for highlighting the implementation of molecular dynamics simulations in the context of plasma sputtering deposition: Zr x Cu 1-x metallic glass and AlCoCrCuFeNi high entropy alloy thin films deposited onto silicon. Effects of depositing atom kinetic energies and atomic composition are studied in order to predict evolution of morphologies and atomic structure of MD grown thin films. Experimental and simulated X-ray diffraction patterns are compared.

Introduction

Plasma sputtering technology is widely used both in research and in industry (Mattox 1989;Safi 2000;Kelly 2000;[START_REF] Musil | Reactive magnetron sputtering of thin films: present status and trends[END_REF]). Nevertheless, this fairly simple technology is continuously being improved [START_REF] Helmersson | Ionized physical vapor deposition (IPVD): A review of technology and applications[END_REF][START_REF] Alami | High power pulsed magnetron sputtering: Fundamentals and applications[END_REF]Sakarinos, 2010). A first interest of magnetron sputtering is to build coatings with improved or complementary properties compared to the pristine bulk materials. Essentially, thin film properties are related on composition and microstructures, which in turn depend on the deposition conditions and thus on sputtering plasma properties. A second specific property of plasma sputtering and thus more specifically magnetron sputtering, is it acts as an atom source, able to interact with an inert and/or reactive buffer gas, which is contributing to the plasma ignition and properties. Such a source is thus providing a sputtered vapour (more or less directed) with well defined composition, kinetic energy and angle distributions. Such "initial conditions" determine the building of the coatings, for which a special initial stage is the growth of supported clusters which may have interests in themselves. Studying plasma sputter growth of supported clusters and thin films remains a challenging task for understanding current properties (film morphology, use properties) and predicting new properties (Safi, 2000;Kelly, 2000;[START_REF] Musil | Reactive magnetron sputtering of thin films: present status and trends[END_REF][START_REF] Brault | Plasma deposition of catalytic thin films: Experiments, applications, molecular modeling[END_REF]). Due to the atomic nature of this deposition process, an insight of the growth processes at the atomic scale is thus of paramount importance. Precisely, molecular Dynamics is a simulation technique that is well suited for describing materials at the atomic scale [START_REF] Alder | Phase transition for a hard sphere system[END_REF][START_REF] Alder | Transport processes in statistical mechanics[END_REF][START_REF] Alder | Studies in Molecular Dynamics. I. General Method[END_REF][START_REF] Allen | Computer Simulations of Liquids[END_REF][START_REF] Haile | Molecular Dynamics Simulations[END_REF][START_REF] Rapaport | The Art of Molecular Dynamics Simulations[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF]) and thus well suited for predicting clusters and thin film compositions and morphologies, as, for example issued from plasma deposition and growth experiments. So, numerous approaches have been proposed for describing plasma-surface phenomena [START_REF] Graves | Molecular Dynamics Simulations of Ion-Surface Interactions with Applications to Plasma Processing[END_REF][START_REF] Graves | Molecular dynamics for low temperature plasma-surface interaction studies[END_REF]Wadley, 2001;Bogaerts 2009Bogaerts , 2010[START_REF] Bogaerts | Modeling of plasma and plasma-surface interactions for medical, environmental and nano applications[END_REF][START_REF] Brault | Plasma deposition of catalytic thin films: Experiments, applications, molecular modeling[END_REF]. Focus on plasma sputter growth using MD is intended for illustrating the power of MD and the capability of predicting morphology, structure and composition in agreement with experimental results (Georgieva, 2010[START_REF] Georgieva | Understanding the Surface Diffusion Processes during Magnetron Sputter-Deposition of Complex Oxide Mg-Al-O Thin Films[END_REF]Xie 2013aXie , 2013b). In the next section, MD simulations are presented in the context of plasma sputter deposition, with emphasis on initial conditions. This will be followed by two original case studies which are highlighting the MD predictive capabilities for complex far from equilibrium systems: Zr x Cu 1-x metallic glass (MG) and AlCoCrCuFeNi high entropy alloy (HEA) thin films deposition and growth onto silicon (100) substrate. Effects of depositing atom kinetic energies and atomic composition will be studied in order to correlate with morphologies and atomic structure.

2. Implementing MD simulations in the context of plasma sputtering deposition.

Classical molecular dynamics

Classical molecular dynamics (MD) is a method intended to numerically solve the Newton equations of motion for a given system of particles, for which the interactions are governed by a model describing the forces between these particles. An advantage of molecular dynamics simulations is that systems can be studied with an atomic resolution at short time-and length-scales, down to femtoseconds and angstroms up to nanoseconds, sometimes to microseconds and a few hundred nanometers [START_REF] Hansen | Modeling of metal thin film growth: Linking angstromscale molecular dynamics results to micron-scale film topographies[END_REF]. A very striking advantage in MD simulations is that no assumption is made about the mechanisms operating in the system. Molecular dynamics simulations consist in tracking the solutions of the classical equations of motion, written as: Where m i , i r , i v and i a are the mass, position, velocity and acceleration of atom i in a defined coordinate system, respectively. The force F i acting on each atom i accelerates the particle i, which in turn leads to a new position, velocity, and acceleration. Forces are derived as the negative gradients of the potential energy (interaction potential) U with respect to atomic positions:

( ) N 3 2 1 ri i r , , r , r , r U F  ∇ - = (2)
The potential energy, can be obtained using either empirical potential energy expressions, semiempirical methods, or ab-initio approaches. In MD simulations, the calculations of the potential energy and force are the most computer time consuming parts. Once the initial conditions and interaction potential are defined, Equation ( 1) and ( 2) can be integrated to follow the time evolution of the atoms in response to the applied forces. The trajectory of all atoms can be thus obtained as solutions of Equation ( 2), in the form of positions and velocities of all the atoms at each time step. The velocity Verlet algorithm [START_REF] Swope | A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters[END_REF] is used for integrating Newton equations. For maintaining the substrate at the desired temperature (T s = 300 K in all simulations here) a Berendsen thermostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF], is applied for a specified duration: the thermal relaxation time. Typical thermal relaxation times range from 0.1 to a 50 ps [START_REF] Hou | Deposition of AuN clusters on Au(111) surfaces: I. Atomic-scale modelling[END_REF] per impact on the surface.

Interaction potentials

In classical MD simulations, the forces between atoms in the system are derived from a potential energy function; this functional form is often based on a quantum mechanical treatment of the system.

The more fundamental quantum mechanical treatment is simplified, by the use of various approximations, the values of which are taken, either from first principle calculations, or from fits of the model to experimental data.

In the present work, Lennard-Jones potential are used to describe the interactions of the Si-(ZrCu) and the Si-(AlCoCrCuFeNi) systems [START_REF] Zhang | Molecular simulation study of nanoscale friction for alkyl monolayers on Si (111)[END_REF][START_REF] Zhen | Calculation of the Lennard-Jones n-m potential energy parameters for metals[END_REF][START_REF] Graves | Molecular dynamics for low temperature plasma-surface interaction studies[END_REF]. While pair potential parameters for the corresponding compound materials are not directly available, the Lorenz-Berthelot mixing rule is used. [START_REF] Allen | Computer Simulations of Liquids[END_REF][START_REF] Graves | Molecular dynamics for low temperature plasma-surface interaction studies[END_REF]. More generally, for any kind of potential we can use, in a first approximation, the following general rules can be established: for length terms we consider arithmetic averaging and for energy terms, geometric averaging. One could guess that more precisely defined potentials were available, which is not the case for all of these metal -silicon substrate interactions. Extracting « true » potentials from ab-initio simulations is out of the scope of this article.

Moreover successful comparison with experiments (Xie, 2013a(Xie, , 2013b) ) allows using this metal substrate interaction potential for understanding effects of incoming atom kinetic energy and composition on morphology and structure of grown thin films.

Metal atom interactions can be well described by the so-called embedded atom model (EAM) [START_REF] Foiles | Contributions of the embedded-atom method to materials science and engineering[END_REF]Maisuradze, 2003). Such a potential, first introduced by Daw, Baskes and Foiles [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][START_REF] Daw | The embedded-atom method: a review of theory and applications[END_REF][START_REF] Foiles | Application of the embedded-atom method to liquid transition metals[END_REF][START_REF] Foiles | Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys[END_REF], is non-pairwise in the sense that it is based on concepts coming from density functional theory which stipulates, in general that the energy of a solid is a unique function of electron density. The embedded-atom method uses the concept of electron (charge) density to describe metallic bonding. Essentially, each atom contributes through a spherical, exponentially-decaying field of electron charge, centred at its nucleus, to the overall charge density of a system. Binding of atoms is modelled as embedding these atoms in this "pool" of charge, where the energy gained by embedding an atom at location r is some function of local density.The EAM parameters of Al, Co, Cr, Cu, Fe, Ni, Zr are given by Zhou et al (Zhou, 2001[START_REF] Zhou | Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers[END_REF] and Lin et al [START_REF] Lin | Computational study of the generation of crystal defects in a bcc metal target irradiated by short laser pulses[END_REF].

Silicon substrate is described using Tersoff semi-empirical three-body potential (Tersoff 1988a(Tersoff , 1988b[START_REF] Tersoff | Modeling solid-state chemistry: Interatomic potentials for multicomponent systems[END_REF]. This potential has been successfully used to investigate the structural, thermal vibration and surface properties of Si.

The EAM and Tersoff parameters are also implemented in the MD LAMMPS software [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF][START_REF] Plimpton | Computational aspects of many-body potentials[END_REF].

There exist many other semi-empirical analytical potentials. Also suited to metals are tight binding potential (Rosato,1989;[START_REF] Sutton | Long-range Finnis-Sinclair potentials[END_REF][START_REF] Cleri | Tight-binding potentials for transition metals and alloys[END_REF][START_REF] Treglia | Alloy surfaces: segregation, reconstruction and phase transitions[END_REF][START_REF] Goyhenex | Theoretical insight in the energetics of Co adsorption on a reconstructed Au(111) substrate[END_REF]. For covalent systems, reactive empirical bond order potentials -REBO and Adaptive Intermolecular REBO - [START_REF] Brenner | Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films[END_REF][START_REF] Brenner | The art and science of an analytical potential[END_REF][START_REF] Stuart | A reactive potential for hydrocarbons with intermolecular interactions[END_REF][START_REF] Ni | Compression of Carbon Nanotubes Filled with C60, CH4, or Ne: Predictions from Molecular Dynamics Simulations[END_REF][START_REF] Nikitin | Hydrogenation of Single-Walled Carbon Nanotubes[END_REF], are now widely used. A new reactive semi empirical ReaxFF potential is growingly attracting attention [START_REF] Van Duin | ReaxFF: A Reactive Force Field for Hydrocarbons[END_REF][START_REF] Nielson | Development of the ReaxFF Reactive Force Field for Describing Transition Metal Catalyzed Reactions, with Application to the Initial Stages of the Catalytic Formation of Carbon Nanotubes[END_REF]Neyts, 2012).

Initial conditions for releasing sputtered atoms.

Plasma sputtering deposition consists in producing an atom source sputtered from biased solid targets by the created plasma ions, typically argon. The sputtered atoms (or molecules, clusters) will lose energy through collisions with the gas phase along the pathway to the substrate. We use a simple model of particle slowing down along straight-line trajectories, subject to a continuous energy loss process, which allows calculating the energy distributions of the particle flux as a function of the distance to the source. The mean kinetic energy of sputtered atoms is calculated according to the modification of Thompson formula, as the sputtered atoms are travelling through a gaseous medium [START_REF] Meyer | Thermalization of sputtered atoms[END_REF]:

3 coh 2 2 1 Ar coh E E 1 E E E E 1 ) E ( f       +         γ + - ∝ + (3) ( ) ( ) [ ] g B i f g B F T k E / E ln n exp T k E E + - = (4) E f /E i = 1-γ/2
is the ratio of energies after and before a collision, [START_REF] Gras-Marti | Slowing down and thermalization of sputtered particle fluxes: Energy distributions[END_REF] with

2 s g s g ) m m ( m m 4 + = γ
where m g and m s stand for plasma forming atom (argon here) and sputtered atom masses, respectively, and n=dpτ/k B T g . E is the energy of the sputtered particles as they leave the target, T g is the sputtering gas temperature, n is the number of collisions that take place in the gas, d is the travelled distance, p is the sputtering gas pressure, and τ is the collision cross section assuming hard core interactions. To calculate the energy loss (E f ) of sputtered atoms with the gas atoms, a Maxwell-Boltzmann (MB) distribution at T g is fixed for the gas (in our conditions T g =300 K, Eg = k B T g ). Because we search for the complete energy distribution of sputtered atoms, for each E g in the MB gas distribution, the energy loss is calculated for a fixed value of the kinetic energy E of a sputtered atom. This is repeated for each E in the Thompson distribution and weighted by the collision probability, which is simply the convolution of f(E) and the MB distribution at T g . [START_REF] Meyer | Thermalization of sputtered atoms[END_REF]. So, initial velocities or equivalently kinetic energies can be obtained from such a distribution using a rejection method [START_REF] Press | Numerical recipes[END_REF]. The velocities of the sputtered atoms can also be sampled from a Maxwell-Gaussian distribution with a mean kinetic energy corresponding to the energy loss calculated with Equation 4 with randomly selected incident angles. Finally Monte-Carlo solvers as SIMTRA [START_REF] Van Aeken | The metal flux from a rotating cylindrical magnetron : a Monte Carlo simulation[END_REF][START_REF] Depla | Magnetron sputter deposition as visualized by Monte Carlo modeling[END_REF] are able to calculate the energy and angle distribution function in three dimensions in the real plasma sputtering system. When the distribution function is found, a rejection method can be used to sample the kinetic energies and incidence angles of the impinging particles on the substrate.

When initial velocities and incidence angle are known, initial positions in the horizontal plane are randomly sampled in uniform distribution at a vertical position a bit greater than the atom -surface and atom-atom cut-off distances. The initial vertical position can be gradually changed as the film is growing for preventing section of position inside the film. For mimicking deposition the atoms should be released to the surface one after each other. The effect of argon ions is neglected in the present simulations, so no argon ions are released to the surface during deposition.

Limitations and capabilities of MD deposition studies.

The time interval between each released atom should be chosen greater than the energy relaxation of the substrate after the collision, which ranges from 1 to few tens of picoseconds [START_REF] Hou | Deposition of AuN clusters on Au(111) surfaces: I. Atomic-scale modelling[END_REF][START_REF] Brault | Platinum nanocluster growth on vertically aligned carbon nanofiber arrays: sputtering experiments and molecular dynamics simulations[END_REF]. It is necessary to be aware that atom fluxes in MD are unphysically large when comparing with experiment. When considering MD simulations, sputtered atoms are released here every τ=2 ps on a 10 x 10 nm 2 surface, this leads to a flux of 10 24 at. cm -2 s -1 . We have checked that increasing this time does not change the results. This should be compared to a typical flux encountered in sputtering deposition: say, 10 15 at.cm -2 s -1 (around 1µm.h -1 ). The first and important consequence is there is no possibility for predicting a deposition rate. Nevertheless, when recalling that excess energy is correctly described using an appropriate thermostat, other information as sticking coefficient, structure, and morphology are reachable. Moreover information as radial density functions (RDF) and simulated xray diffraction (XRD) patterns can be readily compared to experiments [START_REF] Hansen | Modeling of metal thin film growth: Linking angstromscale molecular dynamics results to micron-scale film topographies[END_REF]Xie, 2013aXie, , 2013b)). The second main limitation is the possible long time diffusion occurring in growth process on surface. Some strategies have recently been successfully carried out including Monte-Carlo procedure coupled to MD simulations [START_REF] Blackwell | Modeling evaporation, ion-beam assist, and magnetron sputtering of thin metal films over realistic time scales[END_REF][START_REF] Neyts | Combining molecular dynamics with Monte Carlo simulations: implementations and applications[END_REF]. Such investigations are out of the scope of this article, because in the present studies, individual atomic motions are much localised for HEA and metallic glasses. Indeed, we have the chance that the atom size difference acts for preventing large diffusion leading to the (near)-amorphous phase. The used potentials have the effect of normally leading to crystalline phase for individual atoms (see ZrCu with large Cu of Zr content). Moreover in the present simulations, we have carried out simulation on 100 and 17 Å substrate widths without any significant differences. We think this comes from the size difference between atoms. For atomsubstrate interaction, the cut-off length is always less than half the box size: this prevents the size effects when applying periodic boundary conditions. More precise description of the metal-silicon will require ab-initio determination of the interaction potentials.

Case study 1: Zr x Cu 100-x (10 ≤ x ≤ 90) thin films deposited on a silicon (100) substrate: effect of the deposited atom kinetic energy

Recently, there has been a huge interest in the atomic-level structure and structure-property relationship in metallic glasses (MGs). These materials have been studied for 40 years because of their promising properties belonging to both metals (electron, heat conductivity, ductility…) and glasses (hardness…) [START_REF] Li | Formation of Bulk Metallic Glasses and Their Composites[END_REF][START_REF] Masumoto | Recent progress of amorphous metallic materials[END_REF]. To stabilize an amorphous phase in metallic alloys, atomic diffusion must be prevented. This could be achieved by setting up with the chemical composition (mixing of elements with different atomic sizes) or by freezing a low ordered phase during the synthesis process [START_REF] Kawamura | Fabrication of bulk amorphous alloys by powder consolidation[END_REF]. It has been shown that deposition of thin films by condensation onto cold substrates allows stabilizing low ordered structure in such systems.

As an example, ZrCu alloys have attracted interest in recent years, due to their bulk metallic glass properties [START_REF] Cheng | Atomic-level structure and structure-property relationship in metallic glasses[END_REF][START_REF] Wang | Bulk metallic glass formation in the binary Cu-Zr system[END_REF][START_REF] Yu | Superior glass-forming ability of CuZr alloys from minor additions[END_REF][START_REF] Pǎduraru | An interatomic potential for studying CuZr bulk metallic glasses[END_REF][START_REF] Mei-Bo | Binary Cu-Zr bulk metallic glasses[END_REF], and as amorphous alloy films for their mechanical [START_REF] Xu | Bulk metallic glass formation in binary Cu-rich alloy series-Cu 100-x Zr x (x= 34, 36, 38.2, 40 at.%) and mechanical properties of bulk Cu 64 Zr 36 glass[END_REF][START_REF] Musil | Structure and mechanical properties of magnetron sputtered Zr-Ti-Cu-N films[END_REF][START_REF] Das | Work-Hardenable" Ductile Bulk Metallic Glass[END_REF] and superconductivity properties [START_REF] Karpe | Influence of deposition conditions and ion irradiation on thin films of amorphous Cu-Zr superconductors[END_REF]. [START_REF] Dudonis | Synthesis of amorphous Zr-Cu alloys by magnetron cosputtering[END_REF] prepared thin films with Zr composition in the range of (5 ≤ x ≤ 95) by using high working power (490 W and 1380 W on Cu and Zr targets, respectively) during magnetron sputter deposition.

Numerous theoretical studies have also been conducted on ZrCu systems. Sha et al employed atomistic methods (among them MD) for studying ZrCu MGs forming conditions [START_REF] Sha | Statistical composition-structure-property correlation and glass-forming ability based on the full icosahedra in Cu-Zr metallic glasses[END_REF](Sha, , 2012)).

Almyras et al investigated the microstructure of Zr 35 Cu 65 and Zr 65 Cu 35 MGs and found that these systems consist of small touching and/or interpenetrating icosahedral-like clusters which results in "supercluster" (SCs) satisfying the system composition [START_REF] Almyras | On the microstructure of the Cu 65 Zr 35 and Cu 35 Zr 65 metallic glasses[END_REF]. They thus claimed that seeking the equilibrium configuration for interpenetrating icosaedron-like clusters allows the prediction of the MG microstructure. While bulk amorphous structure is known to be formed under specific synthesis conditions, the dynamics of growth of ZrCu amorphous thin film has not been so much studied using MD. A better understanding of thin film growth can be thus achieved via MD simulations.

Mean kinetic energies are determined using Equations ( 3), ( 4) and the corresponding averaging procedure summarized in section 2.3. The aim is to reproduce the three different experimental situations summarized in Table 1, for which MD growth runs are simulated. This is intended for determining the role of both the composition and the deposited atom kinetic energy on the formation of the expected amorphous phase. Figure 3 shows the intensity and position of the XRD θ-2θ peaks as a function of the Zr metal content.

Again at the Zr concentration of 90%, Zr x Cu 100-x films are crystalline which agrees with snapshot in figure 1. As the Zr concentration increases from 10% to 80%, a shift of the peaks towards small angles (table 2) is observed. At the same time, all of these patterns consist of broader peaks, indicating a low ordered structure. In Zr 90 Cu 10 , Cu atom with the small diameter can enter the gap between the Zr atoms, which has no or only little effect on the overall structure. While for Zr 10 Cu 90 , even a small content of Zr, due to a larger diameter, disturbs the crystalline structure of Cu. Clearly, the film exhibits two grains: one in orientation marked as "A" and the other in orientation marked as "B". At the beginning of the deposition process, the formed stable nuclei are crystalline with a random orientation distribution. As a result, the growth and coalescence of small clusters leads to their coarsening. In this simulation the energy of deposited atom is not high enough to overcome the diffusion activation energy and restructuration between clusters cannot occur. Therefore the boundary is immobile through all the film thickness. As the depositions goes on, the clusters have a similar growth rate but with different growing direction, as a result, two grains are formed with different lattice orientations.

For the condition 2 (medium-high energies, < E Zr > = 7.65 and < E Cu > = 6.67 eV ), the morphologies of films deposited on the Si(100) are displayed in Figure 5. To precisely study the appearance of the amorphous phase, the range of Zr composition is varied between 3 to 95%. The films grew to approximately 4 to 9 nm when increasing Zr concentration. In this condition, for comparison with set 1 of simulations, there is no grain boundary, and the film is single crystalline for Zr 90 Cu 10 and Zr 95 Cu 5 .

For the highest kinetic energies, <E Zr > = 12.6 eV and <E Cu > = 9.61 eV, the films exhibit a more compact structure than for previous conditions, with a thickness of approximately 4 ~ 7.5 nm.

Corresponding snapshots of the films deposited on the Si (100) are displayed in Figure 8. Based on these simulations, the total RDF is calculated for each case and is presented in Figure 9. The same trends as previously reported are observed when increasing Zr concentrations. In the RDF of Zr 90 Cu 10 , four peaks arise up, evidencing a well crystalline structure (Figure 9). From Zr 10 Cu 90 to Zr 80 Cu 20 , the peaks become broad and split evidencing amorphous structures.

Figure 9: Total RDFs for Zr x Cu 100-x alloys.

Figure 10 shows the intensity and position of the XRD θ-2θ peaks as a function of the Zr metal content. In this simulation, at the Zr concentration of 90%, the XRD peak is sharp showing the well crystalline state, which agrees with the RDF presented in figure 9. As the Zr concentration increases from 10% to 80%, the peak is broadening indicating an amorphous structure and the peak position shifts towards small angles (Table 4). When comparing these calculated XRD patterns with lower depositing energy, it is seen that Zr 90 Cu 10 in both sets show a crystallized structure, while for the other ratios an amorphous structure is formed. This indicates the phase and structure of alloy is not much affected by the metal atom kinetic energy in the "high energy" domain. However, Figure 11 shows a specific orientation of the crystal structure of Zr 90 Cu 10 for simulation sets 2 (A) and 3 (B). It can be seen that snapshot A displays a well crystallized structure, while snapshot B exhibits a crystalline structure capped by an amorphous structure, 1.5nm high on the film top.

During the deposition process, the impinging atoms make impacts with the growing film at the surface.

The kinetic energy of depositing atoms is transferred to the atoms near the impact site and the temperature of the system will locally rise. In order to mimic isothermal growth conditions, Berendsen thermostat was used in an intermediate region above the fixed Si(100) substrate, for energy dissipating and maintaining the substrate at a desired temperature. For set 2, the large kinetic energy of depositing atoms can be properly dissipated, facilitating the atoms hopping to the sites of stronger binding and forming a continuous crystalline structure. While for set 3 with a higher kinetic energy, the excess energy of the incoming atoms is perturbing the lattice and physically destroying the lattice, limiting the reconstruction.

Case study 2: Molecular dynamics simulation of AlCoCrCuFeNi high entropy alloy thin film deposition onto silicon

High Entropy Alloys (HEA) are metallic compounds containing six to 13 elements with a concentration in the range 5 to 35% for each element [START_REF] Yeh | Recent progress in high-entropy alloys[END_REF][START_REF] Zhang | Microstructures and properties of high-entropy alloys[END_REF]. In these alloys, because of the high entropy of mixing, formation of brittle intermetallic phases are avoided and multielement solid solutions are stabilized (fcc and/or bcc) [START_REF] Yang | Prediction of high-entropy stabilized solid-solution in multi-component alloys[END_REF]. The combination of numerous metallic elements with different sizes induces a sluggish atomic diffusion, leading to nanocrystalline or amorphous structures. Bulk HEAs have been extensively studied for their excellent properties and performance, including outstanding strength and resistance to oxidation at high temperatures, anti-adhesion, corrosion resistance, high temperature stability, hydrophobicity, high stiffness, strength and toughness, high hardness and high-strain-rate superplasticity [START_REF] Hsu | Alloying behavior of iron, gold and silver in AlCoCrCuNi-based equimolar high-entropy alloys[END_REF][START_REF] Wu | Adhesive wear behavior of AlxCoCrCuFeNi high-entropy alloys as a function of aluminum content[END_REF][START_REF] Lin | Evolution of microstructure, hardness, and corrosion properties of highentropy Al0.5CoCrFeNi alloy[END_REF] . However, HEAs are quite difficult to synthesize as bulk materials and, for many applications, transferring their exceptional properties to a bare material may be of interest. Few studies on HEA thin films deposited by magnetron sputtering have been published (Yeh, 2004;[START_REF] Chang | Influence of substrate bias, deposition temperature and post-deposition annealing on the structure and properties of multi-principal-component (AlCrMoSiTi) N coatings[END_REF][START_REF] Tsai | Diffusion barrier properties of AlMoNbSiTaTiVZr highentropy alloy layer between copper and silicon[END_REF]Chen, 2004[START_REF] Chen | Nanostructured nitride films of multi-element high-entropy alloys by reactive DC sputtering[END_REF][START_REF] Dolique | Complex structure/composition relationship in thin films of AlCoCrCuFeNi high entropy alloy[END_REF][START_REF] Dolique | Thermal stability of AlCoCrCuFeNi high entropy alloy thin films studied by in-situ XRD analysis[END_REF] starting from alloyed targets formed by melting or casting, or from mosaic ones [START_REF] Dolique | Complex structure/composition relationship in thin films of AlCoCrCuFeNi high entropy alloy[END_REF][START_REF] Dolique | Thermal stability of AlCoCrCuFeNi high entropy alloy thin films studied by in-situ XRD analysis[END_REF]. HEAs can be viewed as atomic-scale composites.

Their final properties come from 1) the combination of properties of their constituting elements and 2) the stabilized solid solutions and structures. Thus for the design of new alloys dedicated to a given application, prediction of structural features of the thin film would be of particular interest. Despite this applicative goal, understanding of the growth mode of such complex multielement films is a challenging research subject. Available literature on metallic alloy thin film growth allows thus carrying out simulation of the process at the atomic level.

The present section is devoted to molecular dynamics simulations for studying growth of complex HEA thin films. The integration time-step is chosen to be 1fs. The HEA atoms impinge on Si(100) substrate with velocities sampled in a Maxwell-Boltzann distribution leading to an incident mean kinetic energy of 1eV. Two substrate sizes are chosen for the subsequent simulations: 1) Si( 100 During MD simulation, 10000 atoms are released to the Si surface for realizing three deposition sets.

A random number x is generated for selecting the nature of the atom to be released, where 0 ≤ x < 100.

Taking Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 as an example, Al is chosen

if 0 ≤ x < 2, Co if 3 ≤ x < 11, Cr if 11 ≤ x <
43 and so on. With this method the desired number of different atom types can be properly selected according to the expected film composition. The outputs of simulations will inform about whether targeted composition is reached. Table 5 gives the number of atoms really living in the simulated film and its composition. The compositions of the simulated samples are close to the targeted ones. 13. In all cases, there is no continuous film growth and the deposits grow as various sized clusters, more or less meandering, due to the limited number of impinging atoms and the large simulation cell size. The cluster mean height is around 3 ~ 3.5 nm for 10000 released atoms. The RDF for the three samples was then calculated and is plotted in Figure 14. The peaks represent the successive neighbour distances of the atoms, which can be compared with the corresponding values in crystalline phases of each element (see Table 6).

Figure 14: Total RDF of the HEA deposit at different compositions.

∑ = >= < n , 1 i i i neighbour R c r
according to the element crystal values presented in Table 6, where c i is the atomic percentage of the i th component and R i is the neighbour distances of i th HEA element in pure crystal. By analyzing the average neighbour distances value of HEA in Table 7, we can see that the first neighbour distance is almost the same in these three samples and around 2.45 Å. The expected 2 nd neighbour distance does not appear in the simulated thin film. This can be explained by the hard ball model [START_REF] Kao | Rapidly solidified structure of alloys with up to eight equalmolar elements-a simulation by molecular dynamics[END_REF] and attributed to the effect of mixing numerous elements with different atomic sizes as depicted in Figure 15. The circles represent the first, second, third, fourth and fifth shells, respectively. Because the system is compact, the first shell will be stable or only change a little bit due to the presence of neighbour atoms with different atomic radius. However the second and third shells are indistinguishable because of the large fluctuation in occupation of different sites. This introduces the disappearance of the second peak in RDFs. 8). As can be seen from figure 16 and 17, peak intensities are very different between simulations and experiments. However, the three experimental peak intensities are not comparable among them, due to different set of XRD measurements. Thus no conclusion between experimental and simulated peak intensity evolution can be drawn. This shows that such comparison can only be qualitative.

4.2 MD Simulations of HEA film growth on a reduced size Si(100) substrate.

In the previous section, we found that HEA clusters are formed after 10000 HEA atoms being released on the substrate. If we increase the number of the deposit atoms, we can expect to evolve towards a continuous film morphology. But this will cost too much calculation time. Another method is to reduce the width of the box: the width of the new substrate cell is thus set to dimensions (17×17) Å 2 . The length in x and y direction is thus 1/6 of the previous size. In other words, if the same height is reached, the volume of the new simulation box is 1/36 of the previous size box.

The interaction cut-off distances are set to 8Å. The width of the system was thus chosen in a way that each dimension would be at least twice the cut-off. So boundary effects are expected to be limited.

With this smaller width, if the velocity of deposited atom is too large along x and y directions, it will pass through the periodic boundary several times. This will affect the whole system and may introduce artefacts in the morphology. Figure 18 displays the trajectory of a typical Cu atom in the process of deposition in sample 1. The small blue balls and the blue line represent the position history of the moving Cu atom. Each position is separated by a 4 ps travel time and Cu atom can be captured by the nearby surface atoms with a fast relaxation. This typical trajectory shows that this small substrate size can be used for the MD simulations of HEA film growth without any effect of the periodic boundary condition. It can also be noticed that such systems (HEA) are expect not to form large range order, due to size differences and thus preventing long range diffusion leading to segregation.

Figure 18: Projected Cu atom during deposition (cross section of the cell along z). Trajectory time elapsed between each position is 4ps.

During MD simulation, 2000 atoms are released for realizing three deposition sets with the same compositions as in previous section. These sets of calculation are equivalent to releasing 70000 atoms on the previous larger surface (100x100 Å 2 ), so it is equivalent to observing the film at a deposition time approximately seven times larger. Table 9 gives the real composition of HEA thin film in the experiments, the number of atoms really living in the simulated films and its composition. The films grew to approximately 10 nm. Snapshots of the films with three different compositions deposited on Si (100) are presented in Figure 19. It can be clearly seen that a crystalline structure is formed in sample 1 and sample 2, above an amorphous HEA layer. While for sample 3, an amorphous structure is formed on the overall thickness.

Actually, throughout the process of atom deposition, there exist a competition between the dense packing (for lowering enthalpy) and the randomness of atoms (for increasing the configurational entropy) driven by thermal energy. In each step during simulation, each jumping atom or in-situ shaking is to seek its minimum-energy position under the rule of minimizing energy. But, at the same time, the kinetic energy or momentum also drives atoms to jump, which is against the tendency towards the minimum energy. With the same simulation conditions, the "engine" of this competition, enthalpy and entropy, is determined by the only varied parameter, the composition of the HEA. In sample 1 and 2, the dense packing becomes dominant in making the atomic configuration more ordered: this occurs for low Al content (2 and 3%). In contrast, sample 3 a disordered or amorphous phase tends to stabilize: this occurs for large Al content (39%). This can be interpreted as the highentropy effect. Figure 21 shows the intensity of the calculated X-ray intensities vs. 2-θ of HEA samples. The number of atom in this set of simulation is low, thus it is hard to get XRD pattern corresponding to the experiment. To circumvent this, we have repeated 25 times the original HEA, according to the periodic boundary condition in x and y direction. A system with around 50000 atoms was built for calculating more precisely the XRD patterns. The peak positions are shown in Table 11. The simulated peaks at position 2θ = 42.85˚, 43.6˚ and 44.3˚ correspond to sample S3, S1 and S2, respectively. The S2 being close to the bcc(110)α structure while the S3 peak is rather consistent with fcc(111)β (Xie, 2013b). In this latter case, because the peak is broader than for S1 and S2, the structure is much more short range order. These values are close to the experiment, with a small shift towards smaller angles, which may be caused by the limited size of the simulation box. Besides, with the same number of the deposited atoms, as already noticed, the XRD peak of sample 3 is broader than the others. This is consistent to the disordered structure shown in Figure 19. On the contrary, sample 1, 2 have thin and narrow peaks and the well crystallized structure arises with thickness around 4 nm and 5 nm respectively in Figure 19. The width of the simulated XRD peak is proportional to the "amorphous" or crystallized phase in the simulation. Nevertheless, in experiments (Figure 16) S2 displays the broader peak, while in the present simulations it is S3. This discrepancy is not well explained. Due to limited number of atoms in the simulation compared to experiment, the comparison can only be qualitative. One has to keep in mind that the other atoms in the expanded film are just the duplication of the same unit cell. This is the limitation of this model. Nevertheless, we can use this procedure to get some information about the crystalline or amorphous nature of the films, depending on the composition of the HEA and to compare with the experiment. What is interesting in using small cells is thus to study the structure for larger thicknesses.

Conclusion

Maturity of semi-empirical analytic interaction potential determination allows extensive (parallel) implementation of classical MD simulations for handling complex systems among them thin films growth phenomena. Plasma sputtering deposition being a widely use technology, predicting structures, morphologies and properties at atomic scales is challenging for factory of the future, where integrating computational materials science will help in saving costly experimental tests. Especially the initial steps of growth driving the final microstructure are reachable by MD simulation provided the deposition conditions can be translated as inputs for MD simulations. The energetic of the deposition, i.e. how to correctly describe the energy dissipation through the substrate, can be handled with adhoc velocity rescaling, provided relaxation times are known. The definitive unreachable information is the deposition rate due to the unphysical large incoming particle fluxes compared to sputtering experiments. The present work and the numerous others previously published or being published, constantly demonstrate that the compositions, structure, morphologies and sometimes thin film properties are correctly described using MD simulations.

Cases studies proposed here should be seen as highlighting how MD is able to predict or to be compared to experiments. First, Zr x Cu 100-x thin film growth is studied by MD simulations using initial conditions consistent with magnetron co-sputtering process. The crystallinity of the films was analyzed by calculated RDFs and XRD patterns. The results show that films containing large amount of one of the element (higher than 80%) are crystallized, whereas for intermediate compositions low ordered phase is evidenced. This trend is due to chemical disorder, since both atoms seem to be incorporated in the same lattice (solid solution) which induces distortion of the lattice parameter. Results also show that the morphology of the films is dependent on the mean kinetic energy of the incoming atoms. With higher depositing energy, a compact and uniform film is likely to form. When the kinetic energy of the atoms are above a threshold a top amorphous zone is growing, whereas at lower energy, the morphology is in the form of clusters or film depending on the atomic composition. The more crystalline Zr 90 Cu 10 evolves from a polycrystalline state at low impinging kinetic energy to a pure crystalline one at larger.

Second, growth of complex systems such as HEA clusters or thin films can be studied by MD simulations. Due to short diffusion range, thin film growth has been studied on rather small simulation box while cluster growth was studied on larger simulation box. Thin film structures are shown to be sensitive to the HEA composition. Amorphous films with top crystalline layer are growing for low Al content (2 and 3%). Pure amorphous thin film is observed for the composition with the largest Al content (39%). Clusters remain always amorphous as RDFs display broad peaks. Tentative comparison with experimental XRD patterns is limited to the agreement of the peak positions. Unfortunately no comparison with the peak intensity is possible due to a too high experimental thickness. Certainly comparison between MD simulations and ultra-thin deposited layer would be more fruitful, provided XRD θ -2θ measurements can be possible with a very low film thickness (around 10nm).

  S = 7 cm, P = 3 Pa, Target bias (Zr) = 460 V, Target bias (S = 7 cm, P = 0.03 Pa, Target bias (Zr) = 460 V, Target bias (Cu) = 490 V MD simulation was carried out in a three dimensional cell, which was periodic only along x and y directions. The dimensions of the Silicon (100) substrate are (25×25×10) Å 3 . The first two bottom layers of the substrate are fixed, while the other layers are temperature -controlled layers using a Berendsen thermostat. The Zr and Cu atoms with the different ratios are randomly deposited on the substrate (one atom is released every 2ps). The integration time-step is chosen to be 1fs. The total number of released atoms is 10000, leading to a total simulation time of 20 ns. Each atom is placed in the vacuum slab at a random position from 5 to 7Å above the surface. Zr with composition 3 < x < 95 is selected if a generated random number ξ is less than x/100, else Cu is chosen . The composition of the films is nearly the vapour composition and sticking coefficients for all compositions are close to 90%. The initial velocities are randomly sampled in a Maxwell-Boltzmann distribution with the defined mean kinetic energy.In the first set of simulations, the mean kinetic energies of incoming atoms are low: <E Zr > = 0.13eV and <E Cu > = 0.34 eV. Snapshots of the films deposited on the Si(100) are displayed in Figure1. It is interesting to see that, as the Zr concentration increases, the morphology goes through three different stages: film -cluster -film. As we can see from figure 1, Zr 10 Cu 90 is deposited as a film with a thickness around 4 nm. Then from Zr 20 Cu 80 to Zr 50 Cu 50 , the deposited atoms grow as clusters with height ranging from 6.5 nm ~ 7.5 nm. Finally, for Zr 60 Cu 40 to Zr 90 Cu 10 , a film is formed with thickness ranging from 5 nm ~ 7.5 nm.

Figure 1 :

 1 Figure 1: Snapshots of Zr-Cu coatings deposited on Si(100) substrate at different ratios of Zr and Cu

Figure 2 :

 2 Figure 2: Total RDFs for Zr x Cu 100-x alloys, individual curves corresponding to different alloy compositions are displaced vertically for clarity.

Figure 3 :

 3 Figure 3: Calculated X-Ray intensity vs. 2θ of Zr x Cu 100-x film with different compositions.

Figure 4 :

 4 Figure 4: Snapshot of sputtered Zr 90 Cu 10 with the lowest energies (set 1), blue outline show its original cell in MD simulation and solid line show the crystal lattice orientation. The atom size is reduced for

Figure 5 :

 5 Figure 5: Snapshots of Zr-Cu coatings deposited on Si(100) substrate at different Zr/Cu ratios (the

Figure 6 :

 6 Figure 6: Total RDFs for amorphous Zr x Cu 100-x alloys, 3 < x < 95.

Figure 7 :

 7 Figure 7: X-Ray intensity vs. 2θ of Zr x Cu 100-x film with different compositions.

Figure 8 :

 8 Figure 8: Snapshots of Zr-Cu coatings deposited on Si(100) substrate at different ratios of Zr and Cu atoms (the numbers in brackets are the Zr and Cu atom numbers respectively). •Zr, •Cu, •Si. Mean kinetic energies are <E Zr > = 12.6 eV and <E Cu > = 9.61 eV

Figure 10 :

 10 Figure 10: Calculated X-Ray intensity vs. 2θ of Zr x Cu 100-x film with different compositions.

Figure 11 :

 11 Figure 11: Snapshots of Zr 90 Cu 10 with condition 2 (A) and 3 (B) grown on Si(100) substrate with their specific orientation showing the crystal structure. Solid line in (B) is used to separate the amorphous and crystalline zones.

  ) with dimensions (100×100×20) Å, 2) Si(100) with dimensions (17×17×20) Å 2 . The deposition model of the six atoms is shown in Figure12.

Figure 12 :

 12 Figure 12: Schematic picture of the initial HEA deposition model. •Al, •Co, •Cr, •Cu, •Fe, •Ni, •Si

Figure 13 :

 13 Figure 13: Snapshots of HEA films deposited on Si(100) substrate with large substrate •Al, •Co, •Cr, •Cu, •Fe, •Ni, •Si

Figure 15 :

 15 Figure 15: Hard ball model for HEA (reprinted from[START_REF] Kao | Rapidly solidified structure of alloys with up to eight equalmolar elements-a simulation by molecular dynamics[END_REF]

Figure 16 :

 16 Figure 16: Experimental XRD patterns of HEA film sputtered at 300 K with different compositions

Figure 17 :

 17 Figure 17: The simulated XRD patterns at different compositions.

Figure 19 :

 19 Figure 19: Snapshots of HEA films deposited on Si(100) substrate with small substrate, the dashed line separates the crystalline and amorphous structures. •Al, •Co, •Cr, •Cu, •Fe, •Ni, •Si

Figure 20 :

 20 Figure 20: Total RDF of the HEA deposits at the different compositions of figure 21

Figure 21 :

 21 Figure 21: The simulated XRD patterns at different compositions.

Table 1 :

 1 The calculated mean kinetic energies of incoming Zr and Cu atoms corresponding to three different experimental conditions.

Table 2 :

 2 Position of the main XRD peak for ZrCu alloy in the simulation Unit: 2θ(°)

	Zr-Cu	Zr 10 Cu 90 Zr 20 Cu 80	Zr 30 Cu 70	Zr 40 Cu 60	Zr 50 Cu 50
	main peak	42.45	41.6	40.5	39.7	38.65
	Zr-Cu	Zr 60 Cu 40 Zr 70 Cu 30	Zr 80 Cu 20	Zr 90 Cu 10
	main peak	37.8	37.0	36.25		35.7

Table 3 :

 3 Position of the main XRD peak for ZrCu alloy in the simulation (Unit: 2θ (°)) Cu 97 Zr 10 Cu 90 Zr 20 Cu 80 Zr 30 Cu 70 Zr 40 Cu 60 Zr 50 Cu 50 Cu 40 Zr 70 Cu 30 Zr 80 Cu 20 Zr 90 Cu 10 Zr 95 Cu 5

	Zr-Cu Zr 3 main peak 42.1	42.1	41.7	40.35	39.4	38.4
	Zr-Cu Zr 60 main peak 37.5	36.8	35.95	35.75	35.7

Table 4 :

 4 Position of the main XRD peak for ZrCu alloy in the simulation. Unit: 2θ(°)

	Zr-Cu	Zr 10 Cu 90 Zr 20 Cu 80	Zr 30 Cu 70	Zr 40 Cu 60	Zr 50 Cu 50
	main peak	41.85	41.3	40.2	39.15	38.2
	Zr-Cu	Zr 60 Cu 40 Zr 70 Cu 30	Zr 80 Cu 20	Zr 90 Cu 10
	main peak	37.55	36.75	35.7	35.65

Table 5 :

 5 Targeted and simulated deposited alloy compositions Sample simulation snapshots for Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 , Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 and Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 thin films are reported in figure

	Atoms	Al	Co	Cr	Cu	Fe	Ni
	Sample S1		Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 (8630/10000)
	Deposited number	151	956	2835	3211	1023	454
	Simulation composition	1,7%	11,1%	32,9%	37,2%	11,9%	5,3%
	Targeted composition	2%	9%	32%	39%	11%	6%
	Sample S2		Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 (8466/10000)
	Deposited number	237	2154	1214	1550	1692	1619
	Simulation composition	2,8%	25,4%	14,3%	18,3%	20,0%	19,1%
	Targeted composition	3%	26%	15%	18%	20%	18%
	Sample S3		Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 (8680/10000)
	Deposited number	2977	1084	1222	1688	1146	563
	Simulation composition	34,3%	12,5%	14,1%	19,4%	13,2%	6,5%
	Targeted composition	39%	10%	14%	18%	13%	6%
	4.1 MD simulation of HEA deposition on large Si(100)

Table 6 :

 6 The neighbour distances in single element crystal.

		Lattice	1 st neighbor	2 nd neighbor	3 rd neighbor	4 th neighbor
	Elements
		parameter (Å)	(Å)	(Å)	(Å)	(Å)
	Al (fcc)	a Al = 4.05	2 a Al = 2.86	a Al = 4.05	2 3 a Al = 4.96	2 a Al = 5.73
		a Co = 2.51
	Co (hcp)	c Co = 4.07	a Co = 2.51	2 a Co = 3.09 c Co = 4.07	2 c Co = 5.76
	Cu (fcc)	a Cu = 3.61	2 a Cu = 2.55	a Cu = 3.61	2 3 a Cu = 4.42	2 a Cu = 5.10
	Cr (bcc)	a Cr = 2.88	2 3 a Cr = 2.49	a Cr = 2.88	2 a Cr = 4.07	4 11 a Cr = 4.78
	Fe (bcc)	a Fe = 2.87	2 3 a Fe = 2.49	A Fe = 2.87	2 a Fe = 4.06	4 11 a Fe = 4.76

Table 7 :

 7 Average neighbour distances in the HEA.As displayed in Table7, the average values are calculated by

	Elements	1 st neighbor	2 nd neighbor	3 rd neighbor	4 th neighbor
	Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 (average)	2.52	3.19	4.32	4.99
	Sample1	2.44	-	4.27	4.79
	Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 (average)	2.52	3.19	4.19	5.14
	Sample2	2.45	-	4.30	4.83
	Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 (average)	2.64	3.49	4.48	5.3
	Sample3	2.46	-	4.35	-

Table 8 :

 8 2θ positions of first XRD peaks for HEA in the experiment and simulation HEA Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 HEA Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18

	Experiment

Table 9 :

 9 Simulated deposited alloy composition in the small size cell

	Atoms	Al	Co	Cr	Cu	Fe	Ni
	Sample 1		Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 (1660/2000)
	Deposited number	37	165	523	643	189	103
	Simulation composition	2,2%	9,9%	31,5%	38,7%	11,4%	6,2%
	Targeted composition	2%	9%	32%	39%	11%	6%
	Sample 2		Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 (1670/2000)
	Deposited number	59	456	250	282	313	310

Table 10 :

 10 Average neighbour distances in the HEA compared to element crystal values.

		1 st neighbor	2 nd neighbor	3 rd neighbor	4 th neighbor
	Elements
		(Å)	(Å)	(Å)	(Å)
	Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 (average)	2.52	3.19	4.32	4.99
	Sample1	2.48	-	4.08	4.83
	Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 (average)	2.52	3.19	4.19	5.14
	Sample2	2.49	-	4.15	4.88
	Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 (average)	2.64	3.49	4.48	5.3
	Sample3	2.52	-	4.42	-

Table 11 :

 11 HEA Al 2 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6 Co 9 Cr 32 Cu 39 Fe 12 Ni 6 Al 3 Co 26 Cr 15 Cu 18 Fe 20 Ni 18 Al 39 Co 10 Cr 14 Cu 18 Fe 13 Ni 6

	main peak	43,81	44,41	43,77
	Simulation			Unit: 2θ(°)
	HEA Al 2 main peak	43.6	44.3	42.85

2θ position of main XRD peak for HEA alloy in the experiment and simulation

Experiment

Unit: 2θ(°)

Acknowledgements

China Scolarship Council is gratefully acknowledged for grant #2009 60212.


Giovanni Neglia
email: giovanni.neglia@inria.fr

Xiuhui Ye
email: yexiuhui@gmail.com

Maksym Gabielkov

Arnaud Legout

How to Network in Online Social Networks

In this paper, we consider how to maximize users' influence in Online Social Networks (OSNs) by exploiting social relationships only. Our first contribution is to extend to OSNs the model of Kempe et al.

[1] on the propagation of information in a social network and to show that a greedy algorithm is a good approximation of the optimal algorithm that is NPhard. However, the greedy algorithm requires global knowledge, which is hardly practical. Our second contribution is to show on simulations on the full Twitter social graph that simple and practical strategies perform close to the greedy algorithm.

I. INTRODUCTION

The first motivation of any social network is to foster information propagation using social relationships among users. Therefore, it is important to understand how a user can reach a large population, that is how to best exploit users' influence. Domingos et al. [START_REF] Domingos | Mining the network value of customers[END_REF] introduced first the general problem to select individuals to spread information taking into account their influence. Then Kempe et al. [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF], [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF] defined a general optimization frameworkto solve this problem. However, the main implicit assumption made by Kempe et al. is that once the influential individuals have been identified, they can be recruited in order to spread the information of interest. For instance, recruitment can be made with a monetary transaction. Kempe et al. make two main assumptions for the recruitment process: the recruitment budget is limited, and when users are recruited, they cannot refuse (so the recruitment decision is the one of the user who wants to spread information). However, the recruitment process in Online Social Networks (OSNs) is much different, because it is made using friend requests. Such a recruitment has three main specificities that are not covered by the model of Kempe et al.: i) there is no guarantee that the friend request will be accepted; ii) a friend request is cheap and does not require monetary transaction, so dynamic strategies are possible; iii) the metric to quantify influence in OSNs, unlike the one used by Kempe et al., is not only the number of users that relay the information, but also the number of users that receive the information. Cha et al. [START_REF] Cha | Measuring user influence in twitter: The million follower fallacy[END_REF] give a more general discussion about possible metrics of influence in Twitter.

In this paper, we make the following contributions. i) We extend the model of Kempe et al. with the three specificities of OSNs and show that the initial result of Kempe et al. still holds, that is the greedy algorithm is a (1 -1/e) approximation of the optimal algorithm that is NP-hard.

ii) Using the complete social graph of Twitter crawled in July 2012 [START_REF] Gabielkov | Studying Social Networks at Scale: Macroscopic Anatomy of the Twitter Social Graph[END_REF] and consisting of 505 million nodes and 23 billion edges, we show that, if only the degree and reciprocation probability of each node i (respectively d i and r i ) are known and retweet probabilities are homogeneous, the simple strategy to select the nodes with the largest product r i d i performs at most 2.5% worse than the greedy algorithm. Moreover, selecting users at random achieves similar performance when the replication probability of the cascade process is as large as 1% and only requires to know the users IDs. Some similar results were observed by Chen et al. [START_REF] Chen | Efficient influence maximization in social networks[END_REF] and Habiba et al. [START_REF] Habiba | Working for influence: network density and influential individuals[END_REF], but on much smaller graphs.

II. PROBLEM FORMULATION AND ANALYSIS

In this section, we start by summarizing the contribution of Kempe et al. on which we are building [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF], [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF], then we present our extension of the Kempe's model.

A. Kempe's Model

Kempe et al. modeled the propagation of information in a social network using two different discrete-time models. The first one is the order independent cascade model. Nodes with or without the information are respectively called active or inactive. When one node u becomes active at time t, it has one chance to influence (infect) all its non-active neighbors, who may then become active at time t+ 1. From time t+ 1 on, node u is still active but no more contagious. The contagion attempts from new active nodes at t+ 1 are arbitrarily ordered. The probability of success needs to be specified in order to completely describe the model. A quite general case is when u's success probability to infect v depends on the set S of v's neighbors that already attempted to influence v. We denote such probability p v (u, S). In a decreasing cascade model this probability is non-increasing in S, that is p v (u, S) ≥ p v (u, T ) whenever S ⊆ T . This corresponds to the fact that the more nodes have already tried in vain to infect v, the less likely v is to be influenced by other attempts. Starting from an initial set of active nodes A, the process will stop in at most n -1 steps. The main performance metric of interest is the final set φ(A) of active nodes or better its expected size E [|φ(A)|] σ(A). The second model is the general threshold model. In this case, each node has a monotone activation function f v : 2 V → [0, 1] and a threshold θ v chosen independently and uniformly at random from the interval (0, 1]. A node v becomes active at time t + 1 if f v (S) ≥ θ v , where S is the set of active nodes at time t. Interestingly, Kempe et al. [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF] show that the two models are equivalent, in the sense that for any activation functions f v (.), there exist corresponding activation success probabilities p v (.) such that the probability distribution over the final active sets ϕ(A) is the same under both models.

The optimization problem introduced by Kempe et al. [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF] is to choose the initial set A under the constraint that |A| ≤ K so that the expected size of the active nodes' final set is maximized. They show that the problem is NP-hard, but that a natural greedy heuristic reaches a (1 -1/e) approximation factor for the decreasing cascade model (and for the corresponding general threshold model). The greedy heuristic simply incrementally increases the set A starting from an empty set and adding at each time the node v i that maximizes the marginal gain σ(A ∪ {v}) -σ(A). If at each step the selected node is a 1-ǫ approximation of the best node, then the greedy algorithm achieves a (1-1/e-ǫ ′ ) approximation factor, where ǫ ′ depends on ǫ polynomially. The key for proving this result is to show that σ(A) is a non-negative, monotone, submodular function on sets1 , then the conclusion about greedy algorithm's approximation ratio follows from known results on such functions [START_REF] Nemhauser | An analysis of approximations for maximizing submodular set functions-I[END_REF], [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF].

B. Extension of Kempe's Model to OSNs

In the following, we extend the model of Kempe et al. to the specificities of OSNs, and, for the sake of simplicity, we refer to Twitter in our description of the problems. Twitter is one of the largest social networks, but it differs from other social networks, such as Facebook, because it uses exclusively directed edges (arcs) among accounts. Twitter has no notion of bidirectional friendship, but it allows users to follow other users, i.e., to subscribe for their messages. Following does not require any approval from the user being followed. If Alice follows Bob, then Alice is called a follower of Bob and Bob a following of Alice. Twitter users can retweet received tweets, that is forwarding the tweets to their followers. In this paper, we use the notation (V, E) to refer to the Twitter social graph, where V is the set of Twitter users and E is the set of directed edges. We orient the arcs in such a way that they show the tweet propagation direction, e.g., if A follows B, the arc is directed from B to A, because A receives tweets from B. Kempe et al. for their influence maximization problem made the implicit assumption that once the influential individuals have been identified, they can be recruited in order to spread the information of interest. Recruitment is costly (in terms of money or social investment), so the available budget limits the number of individuals to be selected to K.

The recruitment process in OSNs is different from the one described by Kempe. Indeed, a user v is recruited by u when v receives information from u, i.e., in Twitter terminology, v follows u. This recruitment process in OSNs leads to three specificities. First, even if the user follows the most influential individuals, there is no guarantee that they will follow him back (that is, be recruited), we call this specificity the follow back problem. Second, following is a quite cheap operation in OSNs so that more aggressive dynamic strategies are feasible, we call this specificity the dynamic problem. Third, we quantify the influence not only by the number of individuals who actively replicate the information but also who can see the information because they follow the original tweeter or one of the retweeters, we call this specificity the reader problem. In the following, we extend the model of Kempe et al. to tackle these three specificities. a) Follow back problem: We define the influence of a user u as the expected number of users who retweet u's tweets.

We stress that, in our model, a user can be retweeted only if some other users follow him and decide to retweet his tweets. In the following, we show how we can extend Kempe's model to this scenario. Let a node be active if it has read the tweet and decided to retweet it. Then the original influence maximization problem can be rephrased as follows: how to choose K nodes that should initially tweet the message in order to maximize the expected number of retweets. In this case, p v (u, S) is the probability that node v reads and decides to retweet the message tweeted or retweeted from u, given that the nodes in S have already tweeted or retweeted it.

In the original problem formulated by Kempe et al., it is not specified how the K initial users should be infected, i.e., in our language, how they should be convinced to tweet the message. In this paper, we focus on a specific user u 0 that is trying to maximize his influence and cannot reach other Twitter users through some external communication network. Then, u 0 can only carefully select a given set of users to follow and hope that these users will follow u 0 back and will eventually retweet u 0 's tweets. The strategic choice of u 0 is then the selection of his set of followings in order to maximize his influence. We consider for the moment that u 0 makes this choice once and for all at his registration. We observe that Twitter puts a cap to the maximum number of initial followings that is K = 2000 (this limit increases linearly with the number of followers). More formally, let B denote the set of u 0 's followings and let ϕ(B) be the set of nodes that retweet a tweet originally emitted from u 0 . We can write ϕ(B) = v∈V X v , where X v is a Bernoulli random variable, that is equal to 1 iff node v is active at the end of the cascade. Our problem can be formally stated as follows:

Argmax

B E [|ϕ(B)|] subject to |B| ≤ K. (1)
In the same spirit of Kempe et al. [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF], [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF] we assume to know 1) the probability r u that a given user u reciprocates u 0 if u 0 follows u and 2) the probability p v (u, S) that node v reads and retweets the message tweeted or retweeted by u, given that the nodes in S have already tweeted or retweeted it. The knowledge of p v (u, S) is also required for u = u 0 or u 0 ∈ S.

The greedy algorithm for the follow back problem (1) corresponds to the following behavior: user u 0 selects K followings one after the other, maximizing at each step the marginal increment of the function E [|ϕ(.)|]. Our first theoretical result is the following.

Proposition 1. The greedy algorithm is a (1 -1/e) approximation algorithm for follow back problem (1).

Proof: Let (V, E) be the social network's graph without node u 0 . Consider a new graph ( V , Ê), where for each node u in V we add a new node u ′ and a link oriented from u ′ to u. Let V ′ be the set of these newly added nodes and h : V → V ′ the function such that h(u) = u ′ . On this graph we define success probability functions as follows: cascade model on the graph ( V , Ê). u 0 's choice of the set of his followings B ⊆ V corresponds to the choice of the set A = h(B) of initial active nodes in V ′ . Moreover, the probabilities have been defined in such a way that it is possible to couple the two processes so that φ(h(B)) = ϕ(B) + K, where adding K corresponds to the fact that the initial set of active nodes is counted by φ(.) 1) is equivalent to solve the influence maximization problem on ( V , Ê) with the additional constraint that the nodes can only be selected in V ′ . This does not change the property of the function σ(.), that is non negative, monotone, and sub-modular, then the results in [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF], [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF] still hold. In particular the greedy algorithm is a (1-1/e) approximation algorithm for the influence maximization problem defined on ( V , Ê) and then for the problem (1). 2b) Dynamic problem: We now consider a variation of the follow back problem [START_REF] Kempe | Maximizing the spread of influence through a social network[END_REF] where node u 0 is not required to select all the K followings at once, but u 0 can apply more complex dynamic strategies. For example, node u 0 can stop following nodes that do not reciprocate by a given time T and start following new users. In this way, u 0 can follow during a given time window more than K users (but at most K at the same time) and reach in general a larger number of followers (the number can approach K if there are at least K nodes in the network willing to reciprocate u 0 ). This improvement in comparison to the original problem is obtained at the price of a longer time required to select the best followings. The best possible result achievable by u 0 is obtained if we assume u 0 to know a priori which nodes would reciprocate. For each node v, let R v be the Bernoulli random variable indicating if node v reciprocates node u 0 by time T after u 0 starts following v.

pv (t, S) = p v (t, S) for t = v ′ and v ′ / ∈ S, pv (v ′ , S) = r v p v (u 0 , S) and pv (t, S) = r v p v (t, (S ∪ {u 0 }) -{v ′ }) if v ′ = h(v) ∈ S.
w u ′′ 1 = 1 u ′′ 2 w u ′′ 2 = 1 u ′′ 3 w u ′′ 3 = 1 u ′′ 4 w u ′′ 4 = 1 1 1 1 1 1 1 1 1 (b) readers
(A ⊆ φ(A)). It follows that σ(h(B)) = E [ϕ(B)] + K and the problem (
Clearly it holds r v = E [R v ].
We introduce then the following ideal optimization problem:

Argmax B E [|ϕ(B)|] subject to |B| ≤ K and R v = 1 ∀v ∈ B. (2)
The greedy algorithm for this problem at each step selects the reciprocating node that maximizes the marginal improvement of ϕ(B).

Proposition 2. The greedy algorithm is a (1 -1/e) approximation algorithm for the dynamic problem [START_REF] Domingos | Mining the network value of customers[END_REF].

Proof: The proof is analogous to that of Proposition 1. In this case an additional node u ′ is added only for each reciprocating node u, i.e., for each u ∈ V such that R u = 1, and the probabilities can be updated as follows: pv (t, S) = p v (t, S) for t = v ′ and v ′ / ∈ S , pv (v ′ , S) = p v (u 0 , S) and pv (t, (S ∪ {u 0 }) -{v ′ }) = p v (t, S) if v ′ ∈ S, where only the nodes v ′ such that R v = 1 need to be considered.

While the problem (2) requires to know a priori which users are willing to reciprocate u 0 , the greedy algorithm can be implemented online without such knowledge. This practical greedy algorithm operates in steps, where each step has a duration at most equal to T time units. At each step the user follows the node v that brings the largest marginal increase in comparison to the already selected nodes assuming that v reciprocates. If node v reciprocates by time T , node u 0 maintains user v in his list of followings, otherwise u 0 removes v. The algorithm stops when K users reciprocate or when there are no more users to select in the network. It is easy to check that the practical greedy algorithm selects exactly the same users that the greedy algorithm with a priori knowledge of the reciprocating nodes would, but it requires in general a longer time to execute. The reasoning above leads us to conclude that: Proposition 3. The greedy algorithm for the problem formalized in Eq. 2 can be implemented without a priori knowledge of which users reciprocate, and its expected number of retweets is at least (1 -1/e) of the value obtained by any online algorithm where each node can be selected at most once and reciprocation delays of at most T time units are tolerated. c) Reader problem: Now, we define the influence of a user u as the average number of users who read u's tweets, because they follow u or because they follow someone who has retweeted one of u's tweets. This problem can be mapped to a variant of the previous case (where we consider the number of retweets) introducing opportune nodes' weights. We need to change the original graph (V, E) as follows. For each user u ∈ V , we introduce a new node u ′′ and the directed arcs (u, u ′′ ) and (v, u ′′ ), for each node v such that (v, u) ∈ E (see Fig. 1 (b)). We denote V ′′ and E ′′ respectively the set of new nodes and arcs and ( Ṽ , Ẽ) the new graph. By doubling each node, we can separately account for the two roles of a user as a retweeter and as a reader. Going back to the cascade model terminology, at a given time step if node u is active, the corresponding user has retweeted the tweet, and if node u ′′ is active the corresponding user has read the tweet. In order to correctly model the process, we introduce activation success probabilities as follows: pv (t, S) = p v (t, S) for v ∈ V and pv ′′ (t, S) = 1 for v ′′ ∈ V ′′ . We also introduce nodes' weights w v = 0 for v ∈ V w v ′′ = 1 for v ′′ ∈ V ′′ . Let X v be the Bernoulli random variable that indicates if node v is active when the cascade terminates. The number of users that see the tweet is given by ψ(B) = v∈ Ṽ w v X v where B ⊆ V is the set of followings selected by node u 0 . Two different problems can then be defined depending if the set B has to be selected at the begin or can be changed dynamically, similarly to what is done above. The only difference is the fact that the weighted objective function E [ψ(B)] is considered instead of the unweighted one E [ϕ(B)]. Obviously the function E [ψ(B)] is non-negative and non-decreasing, we can also prove the following result.

Proposition 4. The function E [ψ(B)] is submodular.

Proof: We adapt some results of Kempe et al. [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF] relative to the size of the different sets to the case where we consider a weighted sum of the set elements. We need to prove that ψ(B 1 ∪ {z}) -ψ(B 1 ) ≥ ψ(B 2 ∪ {z}) -ψ(B 2 ) for any z whenever B 1 ⊆ B 2 . Let C = φ(B) be the (random) set of nodes active at the end of the cascade starting from the nodes in B. Imagine now to start a new cascade process on the graph activating node z, but taking into account the fact that all the nodes in C have already tried to infect their neighbors. This new cascade is called the residual cascade process and has success probabilities p

(C) v (u, S) p v (u, S ∪ C).
We denote this new stochastic process as S C (z) and the additional nodes in V \ C made active by it as φ C (z). Kempe et al. [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF]Theorem 3] proved that φ C (z) is distributed as φ(B ∪ {z}) -φ(B). Then it holds:

E [ψ(B ∪ {z}) -ψ(B)] = E   v∈φ φ(B) (z) w v   (3)
Consider C 1 ⊆ C 2 , and the corresponding residual processes S C1 (z) and S C2 (z). If we couple the equivalent general threshold models by selecting the same threshold at each node, Kempe et al. [START_REF] Kempe | Influential nodes in a diffusion model for social networks[END_REF]Lemma 3] showed that pathwise

φ C1 (z) ⊇ φ C2 (z). It follows that v∈φC 1 (z) w v ≥ v∈φC 2 (z) w v and then E   v∈φC 1 (z) w v   ≥ E   v∈φC 2 (z) w v   whenever C 1 ⊆ C 2 .
(4) Let us now consider two cascade processes whose initial activation sets are respectively B 1 and B 2 with B 1 ⊆ B 2 , if we couple them as above, we can similarly show that φ(B 1 ) ⊆ φ(B 2 ), then

P (φ(B 1 ) = C 1 , φ(B 2 ) = C 2 ) = 0 whenever C 1 ⊆ C 2 . (5)
We can now wrap-up our intermediate results.

Let B 1 ⊆ B 2 , then E [ψ(B 1 ∪ {z}) -ψ(B 1 )] = E   v∈φ φ(B 1 ) (z) w v   = C1 E   v∈φC 1 (z) w v   P (φ(B 1 ) = C 1 ) = C1 C2⊇C1 E   v∈φC 1 (z) w v   P (φ(B 1 ) = C 1 , φ(B 2 ) = C 2 ) ≥ C1 C2⊇C1 E   v∈φC 2 (z) w v   P (φ(B 1 ) = C 1 , φ(B 2 ) = C 2 ) = C2 E   v∈φC 2 (z) w v   P (φ(B 2 ) = C 2 ) = E [ψ(B 2 ∪ {z}) -ψ(B 2 )] ,
where we have used Eqs. ( 3), ( 4) and ( 5).

From the general results for non-negative, non-decreasing submodular functions, it follows that the greedy algorithm that selects orderly the u 0 's followings that incrementally maximize E [ψ(.)] guarantees a (1 -1/e) approximation ratio. The same result holds in the dynamic case. Due to lack of space we do not define formally the two problems and the corresponding propositions, but we summarize our conclusions as follows.

Proposition 5. The greedy algorithms for the static and dynamic versions of our problem reach a (1-1/e) approximation ratio also when the objective function is E [ψ(.)], the expected number of users who see the tweet.

III. SIMULATIONS ON TWITTER

We have shown in the previous section that for all three specificities of OSNs, the greedy algorithm is a good approximation of the optimal algorithm. As the optimal algorithm is NP-hard, it is a major improvement. Unfortunately, the greedy algorithm requires to know the topology of a social graph as well as all the functions p v (u, S) for every node v, a requirement that is not feasible for the social graph of the size of Twitter. Moreover, greedy algorithms are computationally expensive, because of the inherent cost of evaluating the expected size of φ(A), that can only be estimated by Monte Carlo simulations of the cascade process on the social graph.

In this section, we show using simulations on the full Twitter social graph that simple and practically feasible strategies perform close to the greedy algorithm. We describe the methodology we used to perform our simulations in Section III-A and we discuss the results in Section III-B.

A. Methodology

For our simulations, we considered the simple case when p v (u, S) = p is constant and evaluated different selection strategies on the complete social graph of Twitter as crawled in July 2012 [START_REF] Gabielkov | Studying Social Networks at Scale: Macroscopic Anatomy of the Twitter Social Graph[END_REF]. We had to solve two main problems.

The first problem to solve with the simulations is to manage the large size of the Twitter social graph. Indeed, this graph consists of 505 million nodes and 23 billion arcs and requires roughly 417GB of storage in the form of edgelist. A naive implementation of the simulation would require to load the graph into memory and then use Monte Carlo simulations of the retweet process in order to estimate the objective functions with high accuracy. The followers of the initial node u 0 would retweet with probability 1, then their followers will retweet with probability p and so on until no new node is retweeting. To reduce the memory required to make this computation, we introduce the concept of pruned graphs. A pruned graph is obtained from the original one by sampling each edge with probability p (and with probability 1 -p the edge is removed from the graph). Computing the set of reachable nodes from u 0 on a pruned graph is equivalent to counting the number of retweeting nodes in a specific sample of the retweet random process, but memory requirement is reduced by a factor p (usually p << 1) at the expense of storage increase, because we need to work on multiple pruned graphs (see discussion below) in order to reach the required accuracy.

However, for large values of p the pruned graphs are still large, so to reduce further the memory requirement, we use a two step compression. First we compute the Strongly Connected Components (SCCs) of the graph. Second, we construct a Directed Acyclic Graph (DAG) by abstracting each SCC as a single node and replacing multiple arcs between the nodes with a single arc. Provided that p is quite large, we will observe big SCCs in the pruned graph, thus we can achieve a big reduction in size using our approach. Then we compute the reachability on the obtained DAG and deduce the reachability of the original pruned graph by taking into account the number of nodes in each SCC and the fact that the nodes belonging to the same SCC have the same reachability. This approach decreases the computation time, as well as memory and storage requirements (because of the more compact DAG representation). Both these expedients were proposed by Chen et al. [START_REF] Chen | Efficient influence maximization in social networks[END_REF], but we have been able to make computations on a graph 4 orders of magnitude larger.

The second problem to solve is to determine how many simulations we should run (that is, how many pruned graph we should compute) in order to achieve a given precision for the estimates of the expected number of retweets. To this purpose, we have used two different models that we can only describe shortly here because of space constraints (details are available in a technical report [START_REF] Neglia | How to Network in Online Social Networks[END_REF]). The first model approximates the cascade process with a branching process, where the probability to have a follower with k followers is kq k / k , where q k is the distribution of the number of followers in the graph and k is the average number of followers (see for example [START_REF] Barrat | Dynamical Processes on Complex Networks[END_REF]Chapter 8] for a justification of such an expression). This model requires that different active nodes have different followers and is good only for small values of p. In particular, it may be accurate only when the branching process dies out with probability 1, because otherwise the model predicts that the expected number of active nodes is infinite while this number is obviously limited by the total number N (= |V |) of nodes in the graph. The branching process extinguishes with probability 1 if p k kq k / k < 1, i.e., p < 2 × 10 -4 in the considered Twitter graph. The second model addresses the case for p > 2 × 10 -4 . In this case the branching process theory predicts that the process can still extinguish with a probability p ext that is a decreasing function of p. The intuition behind the second model is to couple the branching process and the actual cascade model and assume that the cascade will reach almost all the N nodes when the branching process does not extinguish and a negligible number of nodes when it does. In particular, given that we are interested in providing upper-bounds for the variability of the process, we simply consider that the number of active nodes is a random variable that is equal to N with probability 1 -p ext and to 0 with probability p ext . Some further refinements of the model lead to the conclusion that the number of samples needed to achieve a reasonable prevision is below 100 for all the values of p we considered, i.e., p = 10 -4 , 10 -3 , 10 -2 , 10 -1 . This result is quite surprising, given the high variability of the degree distribution {q k } (that is power law) and the even higher variability of the skewed distribution {kq k / k }.

With the approach described above we have been able to perform our simulations on the real Twitter social graph with hundreds of millions of nodes and tens of billions of arcs.

B. Results

In this section, we compare using simulations three selection strategies for the user u 0 to select his followers: i) the greedy strategy described above; ii) the high-degree strategy which consists in picking the nodes orderly according to their number of followers (from the largest number to the smallest); iii) the random strategy where followers are selected uniformly at random from the whole set of users.

We have to evaluate the impact of two main parameters: the retweet probability p and the reciprocation probability r. We start by considering the impact of p on the three selection strategies by assuming r = 1, i.e., u 0 is followed by every node it follows. While unrealistic, this assumption is an easy way to assess the impact of p only. Then in a second set of simulations, we relax this assumption and define r as a function of the number of followings and followers of each node.

In our first set of simulations, we consider r = 1 and that u 0 's followers retweet u 0 's tweets with probability 1, while all the other users retweet with probability p ∈ {10 -4 , 10 -3 , 10 -2 , 10 -1 }. Fig. 2 shows the expected number of retweets versus the initial number K of followers u 0 can choose. The average number of followers in the original graph is k ≈ 45, but the effective density, as defined by Habiba et al. [START_REF] Habiba | Working for influence: network density and influential individuals[END_REF], is k p ≈ 4 * 10 -3 for p = 10 -4 (in Fig. 2 (a)). Then, in this case, most of the nodes are not retweeted by any follower, but due to the skewness of the distribution {q k } (the number of followers can be as high as 24,635,412), there are some hubs in the social graph that have an expected number of retweeters significantly larger than 0. The cascade processes from different followers of u 0 do not overlap much (each pruned graph is almost a forest of smalldepth trees with a multitude of singletons), so that the highdegree strategy performs almost as well as the greedy strategy. Due to this structure, the expected number of retweeters significantly increases as the number of u 0 's followers keeps increasing. The random strategy performs poorly, because with high probability the 200 selected followers will be singletons. When p = 10 -3 (Fig. 2 (b)), the cascade originated from the node with the largest degree is already able to reach about 3 * 10 6 users (roughly 1% of the whole social graph) and both the greedy and high-degree strategy select this node first. The other followers selected from u 0 using these two strategies lead to a minor improvement: even adding 199 more followers, the expected number of retweeters increases by only 33%. The random strategy starts paying off because there are much less singletons in the pruned graph. Further increasing p to 10 -2 , the contribution of the first follower is even larger and the contribution of the others even more marginal, as it is shown in the plot in Fig. 2 (c). In fact, a non-negligible SCC appears in most of the pruned graphs and a careful choice of the first follower allows u 0 to have roughly one tenth of the nodes retweeting his tweets (this follower is not necessarily in the largest SCC, but he can reach it). The other 199 followers provide roughly 5% more retweeters. We observe that the effective density is about 0.4, then more than half of the nodes have 0 out-degree/in-degree in the pruned graphs. The greedy strategy and the high-degree one lead to a difference in the number of retweets lower than 2.5%. Interestingly, the probability to randomly pick a node in the largest SCC is now quite high, so the random strategy performs close to the other two strategies, but with a higher variability as shown by the large confidence intervals in the figure. Moreover, the figure shows how the good follower is very likely to be selected among the first 10-20 nodes. The same reasoning allows to explain also the curves in Fig. 2 (d) for p = 0.1. In this case, greedy and high-degree are almost indistinguishable and random has almost the same performance for K ≥ 20.

In our second set of simulations, we consider r v = min{ #followings v #followersv+100 , 1}, where #followings v and #followers v are respectively the number of followings and followers for user v. The rationale behind is that a user with a lot of followers and a few followings is not likely to reciprocate u 0 . We do not claim that this formula has any particular value, apart from allowing us to simply test the effect of heterogeneous reciprocation probabilities. In this case the highdegree strategy selects the nodes according to their effective degree r i d i . Surprisingly, the results for p = 0.01 are qualitatively unchanged as shown in Fig. 3 (a). So the reciprocation probability does not seem to significantly impact the respective performance of the three considered strategies. We have also compared the different algorithms in terms of the expected number of users who can read the tweet. We see in Fig. 3 (b) that the number of readers is obviously much bigger than the number of retweeters, but there is no significant difference in the relative performance of the three strategies.

IV. CONCLUSIONS

In this paper we have considered a user of a social network who tries to maximize his influence through a careful networking strategy. We have shown how greedy algorithms guarantee a good 1-1/e approximation ratio, but much simpler strategies like selecting users with the largest number of followers or even selecting random users may practically reach the same performance on real online social networks.

Fig. 1 (

 1 a) illustrates the graph transformation when r v = r and p v (u, S) = p. We consider now the order independent

Fig. 1 .

 1 Fig.1. Graph transformations. The original nodes are green. The added nodes are red/blue, added arcs are dashed. We have specified the new success probability functions for the simple case when pv(.) and rv are constant and respectively equal to p and r.

1 Fig. 2 .

 12 Fig. 2. Performance of greedy, high-degree, and random strategy on the Twitter social graph for different values of retweet probabilities p. The reciprocation probability is r = 1. Error bars show the 95% confidence intervals.

Fig. 3 .

 3 Fig. 3. Extensions of the simulations by taking into account the probability r that users will follow back (a), and by looking at the number of users who received the tweet instead of the number of users who retweeted it (b).

A set function f (.) is submodular if f (S ∪ {z})f (S) ≥ f (T ∪ {z}f(T ) whenever S ⊆ T and it is monotone if f (S ∪ {z}) ≥ f (S) for each S and z.

Also the results for the case when the greedy algorithm selects at each step a (1ǫ) approximation of the best node can be extended to our case, but for the sake of conciseness, we only refer to the simpler case.

This research is partially supported by Alcatel Lucent Bell Labs in the framework of the ADR Network Science. The authors would like to thank Alonso Silva (Alcatel Lucent Bell Labs), Paolo Giaccone (Politecnico di Torino) and Damien Saucez (Inria) for the helpful discussions.