Ab-Initio Modeling of Adhesive Behaviors at Material Interfaces

In this chapter, density functional theory (DFT) is employed to identify the essentials of adhesive formation at the Al/Ceramic interfaces. Analyses of electronic structure in the DFT used for atomic bond evolution at interfaces, indicate that, adhesion energy mostly dominated by surface energy, may lead to formation of adhesives, which bind two type bulk-materials strongly together, especially for those polar systems involving cubic oxides, carbides and nitrides. In addition, ab-initio molecular dynamics (AIMD) based upon the DFT is adopted, for example, in the simulation of chemical reaction between two contacting reactive slabs: pure aluminum and iron-oxide. This may provide an insight into the dynamical formation of an adhesive (amorphous Al2O3-texture) occurring between two nascent material surfaces if lubricants are not present or are insufficient. Such a texture may bond onto a hard-roller surface as a protective thin film to resist the elevated temperature.


Introduction
Adhesives have been widely used to improve mechanical properties of materials such as hardening, deformation stability of surfaces or thin films for a variety of technological applications [1][2][3][4][5][6].A significant effort focused on exploring nanomechanisms at interfaces where two solid surfaces contacted.For examples, adhesives were critical to successful manufacture for semiconductor circuits.Each deposited layer must adhere strongly to that beneath it to survive thermal stresses.Similarly, a recent report told that, wear of materials usually resulted in a very high economic loss to every country, for instance, economic losses due to the wear would amount to about 1.0-2.5% of the gross national product (GDP) in the USA [7].Therefore, adhesives may have a tremendous effect on the US industry because it determines the lifetime of many working parts, from components of automotive engines to airplane brakes.A particular example is the aluminum-alloys processing which is also relevant to many other metallic alloys, that is, the Al-alloy products are manufactured through a series of processing steps that begin with casting an aluminum ingot.In order to reduce ingot thickness, it passes through a series of rolling operations under high loads at speeds up to 4500 ft/min.In this case, there will take place adhesives (the melting aluminum bonds to steel surfaces as chemical additives).On the positive side, many machining operations rely upon adhesive wear processes (from oil drilling to polishing a Si-wafer).The processing optimization is critical to commercial success.While on the negative side, an interest in dry machining is rising because oils used in the machining are likely evaporated due to frictional heating, and such a vapor impact on the long-term health of machine operators is increasingly concerned.Therefore, the better understanding of adhesion mechanisms may result in positive impacts on many industrial manufacturing.
In this chapter, we will introduce our previous studies of adhesives on aluminum and other alloy surfaces: (1) Adhesion is of particular interest.Although a work has reported the adhesion character at an interface between two pure materials (aluminum-diamond) [8], we will discuss an application of electronic structure to this field between aluminum and various ceramic compounds, to determine factors (e.g., an adhesion energy) that may qualitatively control the adhesion strengths at interfaces because this discussion may approach more to the reality even if such computed factors often show too large deviations from some continuum-scale (finite element) calculations.In addition, some particular adhesives, which bind two type bulkmaterials strongly at the Al/Ceramic interfaces, are discussed.(2) Another simulation from abinitio molecular dynamics (AIMD) based upon density functional theory (DFT) is present to bond breaking and formation, and adhesive transfer of materials over an interfacial area where two different material slabs (with pure aluminum and iron oxide) are compressed and then stretched out.In this study, contact of two reactive slabs may lead to a strong exothermic reaction to form a dynamical thin film bonding onto a hard-roller surface to prevent the surface from the melting.
In the following discussions, although we mainly focus on aluminum/aluminum-oxide samples, without loss of generality, we believe what we discuss here is also relevant to many other materials in which adhesion may occur.

Density functional theory
Macroscopic materials consist of an enormous huge number (~1023 per cm 3 ) of nuclei and electrons, and these particles interact in a complex fashion.With the development of quantum mechanics in the early 1900s, the law governing those interactions was made clear, that is, the so-called many-particle Schrödinger equation.In order to solve this equation efficiently, the Born-Oppenheimer (BO) approximation [9] was usually invoked to separate the ionic degrees of the freedom from electronic ones.Even so, an exact solution for this many-particle problem was still impossible.In the early 1960s, density functional theory (DFT) was proposed as an alternative to the wavefunction-based description of quantum mechanics.The core concept of the DFT was to describe the state of a system not with the complex all-electron wave function (AE-WF), but with the much simpler charge density.Consequently, the DFT may drastically reduce the degrees of the freedom needed to describe the system, yet in principle maintain an accurate depiction of the physics.In recent years, the DFT has become the method of choice for describing large molecular clusters and solid systems.
DFT was based upon theorems developed by Hohenberg, Kohn and Sham [10,11], which stated that the ground-state properties of a system were functionals of the charge density n(r) alone.This charge density, which minimizes the energy functional under a given external potential: () -the interaction between electrons and nuclei, is the ground-state charge density.So the total energy of a system can be expressed as a functional of the charge density, that is, (1) According to the Kohn-Sham ansatz [12], the kinetic term [()] is replaced with a kinetic energy of non-interacting electron gas: [()]; the electron-electron interaction potential [()] is replaced with a classical Coulombic interaction: [()] -the Hartree energy, plus an additional non-classical term: [()] -the exchange-correlation (XC) energy.Therefore, Eq. ( 1) is modified as follows: (2) The advantage of this modification is that the first two and the last terms in Eq. ( 2) can be dealt with a traditional manner; while the third term the exchange-correlation energy, whose form is in principle unknown, can be simplified through several approximations.Eqs. ( 1) and ( 2) are called the N single-particle Kohn-Sham (KS) equations.
In order to solve the total energy [()] in Eq. ( 2) more accurately, a set of Lagrange multipliers , which correspond to the orthonormality of the N single-particle orbitals , are adopted in an effective Kohn-Sham Hamiltonian operator , so the ground-state energy [()] can be obtained through minimizing the charge density n(r), see [11,12].One advantage of this method is that the number of the N single-particle KS equations increases linearly with the system size and number of electrons.
Even if the Kohn-Sham ansatz is exact in principle for ground-state systems, in practice the approximations are necessary because the exact functional form [()] in Eq. ( 2) is still unknown.Consequently, an accurate, yet efficient expression for the XC energy is essential.Even though some scientists were devoting themselves to extending the DFT to treat excited states in materials [13], most of current works mainly focus on systems at the ground state.

The exchange-correlation functional
Several approximations for the XC energy functionals have been developed.The most commonly used approximation in the DFT calculations was the so-called local density approximation (LDA) [14], which was based upon the so-called uniform electron gas.The LDA was an exact approximation for a uniform electron gas in a system with a slowly varying charge density [15].Actually, the charge density n(r) usually behaved like a rapidly varying function of position in some systems.Thus, the corrections to the XC energy based upon the gradient of the charge density were required to reduce the LDA error.
One famous commonly used extension to the LDA was the generalized-gradient approximation (GGA).An example of a GGA functional (PW91) predominantly used in solid-state calculations was that of Perdew and Wang [16].Another modified expression known as the PBE functional has been suggested by Perdew, Burke and Ernzehof [17].Comparing the GGA to the LDA, the former has demonstrated to tremendously improve the accuracy of calculations on solids and molecular clusters [15].

Plane waves and supercell
The method of Plane Waves (PW) basis-sets (orbitals), which was more suitable for calculations of periodic solids in the reciprocal space than others on extending its accuracy to any desired convergence, has been widely used for analyzing adhesion properties through calculations of various electronic structures.To simplify, they can be handled by sampling some special kpoints in the reduced Brillouin-zone to determine the electronic wave function and the electric potential through the filled electronic bands.Two widely used methods were Chadi-Cohen [18] and Monkhorst-Pack [19] samplings.
In a larger supercell, for the wave vector selection, it is often sufficient to only consider a particular value-point-the Γ point because this point is the original position in k-space, and such a choice would confer significant savings with data storage and computational cost.

Pseudopotential approximation
The Plane Waves (PW) method cannot directly be carried out in the KS theorem since the wave functions of ionic-cores are rapidly varying with many nodes and require enormous PW basissets to reach the calculation accuracy.Since chemical reactions are solely related to valence electrons while ionic-core electrons are approximately "frozen" in their original configurations, in practice, to handle the external potential: (), a scheme of pseudopotential approximation (PP) is provided to just reproduce identical curves of valence wave functions in the place sufficiently far from the ionic-core region (R C ).So in the PP framework, the constructed charge densities outside the ionic-core region (R C ) must be identical to the true ones, which is called "the norm-conserving condition" [20].
To compromise the computational demanding in practice, a successful concept of "the ultrasoft-pseudopotential (US-PP)" was introduced by Vanderbilt [21].Blöchl further expanded the US-PP concept by combining ideas from the PP and the LAPW (linearized augmentedplane-wave) methods in a conceptually elegant framework, he called it "the projector augmented-wave method (PAW)" [22].And then, Kresse established an exact formal relationship between the US-PP and the PAW methods [23].

Molecular dynamics
Molecular dynamics (MD) is a methodology which depicts motions of a many-particle system using classical Newton's equations.It has been over 40 years since the first application of MD method created by Alder and Wainwright to a system with hard spheres [24].The MD method is particularly useful for the study of dynamical properties of materials.More important, this study may help people extract physical insights from the many-particles system.

Equations of motion
In MD simulations, classical trajectories of particles in the real space for a finite size system are generated by solving the Lagrange equations numerically.When choosing the Cartesian coordinates, these equations will become Newton's equations, that is, where ( ) is the total potential energy of a system, m i is the mass of the particle i, P i is the momentum of the particle i, F i is the total inter-particle force acting on the particle i, and " ⋅ " and " ⋅⋅ " are the first and the second orders of the time derivative, respectively.The convenience of solving Eq. ( 3) should depend on some practical algorithms.

Algorithms
A standard way of solving Eq. ( 3) is the finite-difference method: Under a given configuration (positions and velocities of all particles) and other dynamical information in a system at the time t, a numerical integration would determine the new configuration of the system at the later time t+Δt (Δt is a time step).Commonly used methods in MD are the Verlet algorithm [25], the Leapfrog algorithm [26] and the Gear predictor-corrector (GPC) algorithm [27].An ideal algorithm should be simple, run fast and require little memory.It should also permit using a longer time step Δt to hold the whole trajectory of the system as close as possible and preserve conservation laws of momentum and energy.According to these, the Verlet algorithm can be the most widely used one.The leapfrog algorithm is essentially identical to the Verlet algo-rithm.The GPC method is often more accurate as well as more complicated than those two algorithms.

Statistical ensembles
The widely used statistical ensembles in MD simulations are often the microcanonical ensemble (NVE), the canonical ensemble (NVT), the isothermal-isobaric ensemble (NPT) and the grand canonical ensemble (μVT).Thermodynamic variables in above parenthesis for each ensemble are fixed in relevant ensembles [28].In MD simulations, the (NVE) ensemble is the best to be realized since all equations of the motion conserve the total energy E and particle's number N. The constant volume V can be fixed using the free or the periodic boundary conditions; The (NVT) ensemble is very commonly used for practical calculations.Therefore, to control the temperature T in this ensemble, there are a number of ways to simulate it: (i) the simplest one is the velocity scaling method [28], but is crude in MD simulations.(ii) the Nosé-Hoover's method is a more conceptually fundamental, but has an adjustable parameter of the thermal inertia in practice [29].
For the (NPT) ensemble, Anderson provided an "extended system" method for a constant pressure MD [30].In this method, a term of "external piston" was added into the Lagrange functional for adjusting a value of external pressure.In another way, Parrinello and Rahman [31] proposed a scheme of "variable cell constant MD" which allowed the simulation box to change the shape as well as the size [31].This method, so far, is the most general constant pressure method in MD simulations.
The (μVT) is very useful for an open system: In this ensemble, the total particle's number N will not be conserved, see [32] for details.

Ab-initio molecular dynamics
The MD simulation based upon the force calculations through empirical potential [33] is called the classical molecular dynamics simulation, while that based upon the force calculations using density function theory is called the ab-initio molecular dynamics (AIMD) simulation [34].
Comparing these two, the former can approximately treat a large system (up to tens of millions of atoms) for long time steps (nanoseconds); while the latter can more accurately simulate a small system of hundreds of atoms, and is better at describing reactions involving the change of atomic bonding types.
In principle, the ab-initio molecular dynamics (AIMD) is an extension of the classical molecular dynamics (MD), which includes the first principle-derived potentials for the calculation of the inter-particle forces and the total energy of a system, see [34,35].

Born-Oppenheimer molecular dynamics
In the AIMD scheme, the total potential energy ( ) of a system has the same physical meaning as the Kohn-Sham energy ( ) by means of the Born-Oppenheimer (BO) approxi- Adhesives -Applications and Properties mation.Since the Kohn-Sham energy depends only upon nuclear positions and the defined hyper-surface for nuclear motions, the Lagrange functional of the so-called Born-Oppenheimer molecular dynamics (BOMD) has the following formula [35] (4) where are the Lagrange-multipliers.So the minimization of Kohn-Sham energy, E KS , is the constraint to the orthogonal basis-orbitals {φ i,j }.And then, Eq. ( 3) for many-particles can be modified by means of Eq. ( 4) as follows: ( where denotes an effective Kohn-Sham Hamiltonian operator for a single-electron.Based upon Eqs.(5-6), a classical molecular dynamics program can be turned into the BOMD program by replacing the energy and the force routines with the corresponding ones from a quantum chemistry methodology [35].The accuracy of the force calculations used in the BOMD depends linearly on the accuracy of minimizing the Kohn-Sham energy E KS .To realize this, an algorithm of the conjugate gradient (CG) minimization was effectively applied to such calculations, see [34].

Car-Parrinello molecular dynamics
The basic idea of the Car-Parrinello molecular dynamics (CPMD) can be seen as exploring a time-scale separation of "fast-electronic" and "slow-nuclear" motions by transforming it into a dynamical framework of "classical-mechanical and adiabatic energies" [36].In classical mechanics, forces on nuclei are obtained through the derivative of a Lagrange functional with respect to nuclear positions.The CPMD scheme suggested that, a functional derivative with respect to particle orbitals, which were interpreted as classical fields, yielded the forces on those suitable particle orbitals given by the Lagrange functional.Thus, the CPMD scheme postulates the following Lagrange functional (7) Like in classical mechanics, nuclear positions, electronic orbitals and the Car-Parrinello's equations of the motion are modified as the following formulae corresponding to the Newton's equations [35] ( where μ is a "fictitious mass" or "inertia parameter" assigned to the orbital degrees of the freedom.The unit of this "mass parameter" μ is the energy multiplies the time square for reasons of dimensionality.Note that constraints within E KS would lead to "force constraints" in Eq. (8).In practice, these constraints would depend upon both the Kohn-Sham orbitals and the nuclear positions through an overlap matrix of basis-functions-the simulated annealing to search for coefficients of the basis-orbitals that minimize electronic energy.In this terminology, a scheme of "low electronic temperature" or "cold electrons" in Eq. ( 9) is close to its instantaneous minimum energy: min ; , that is, in this case, it is close to the exact Born-Oppenheimer hyper-surface.Therefore, a ground-state wave function optimized for the initial configuration of nuclei would stay close to its ground state during time evolution if the system is kept at a sufficiently low temperature.A SHAKE algorithm is used to impose constraints on the system to ensure that the KS orbitals remain orthonormal [34].In practice, it was found that the Car-Parrinello scheme was unnecessary for electronic configurations to be minimized in a coefficient space at each MD time step even though this manner may give errors in force calculations on nuclei.However, it was rather strange to indicate that, the errors of force calculations on nuclei were canceled by the associated errors in depicting electronic motions.According to this, the Car-Parrinello scheme may tolerate force errors on nuclear motions in the AIMD simulation.
In addition, the last term on the right-hand side of Eqs. ( 5) and ( 8) would vanish if basisorbitals {φ i } are independent of nuclear positions {R N }.Comparing the BOMD to the CPMD schemes, the former can be viewed as an alternative to the latter, which may give very accurate coefficients of basis-orbitals through minimizing electronic energy at each MD time step.

An implemented software: VASP
The Vienna Ab-initio Simulation Package (VASP) is a commercial software under the framework of density functional theory (DFT) code that is developed and distributed by a research group at The University of Vienna, Austria [37].VASP uses a plane wave (PW) basis-sets (orbitals) to expand the single-electron Kohn-Sham wave function and the pseudopotentials for describing electron-ion interactions.Three separate approximations to the exchangecorrelation (XC) energy are employed.The local density approximation (LDA) and the generalized-gradient approximation (GGA) are in two flavors: the PW91 functional and the PBE functional [14][15][16][17].In addition to these functionals, several additional routines have been implemented to allow for estimation of transitional state energies and geometries (NEB, metadynamics [38]), vibrational properties, analyses of various electronic structure aspects, etc.For more details regarding the VASP code, see [37].

Density functional study of adhesion at the metal/ceramic interfaces
A fundamental quantity to influence mechanical properties of a material interface is defined as adhesion energy: W sep (also commonly referred to as a "work of separation"), breaking interfacial bonds and reversibly separating the interface into two free surfaces, but neglecting diffusion and plastic deformation [39][40][41].Experimental observations on this dissipative process indicated that, due to plastic deformation, the W sep was proportional to, but smaller than, the work done to cleave the interface, see [42].Thus formally, the W sep is defined in terms of surface and interface energies of respective bulk materials, or to the difference in total energy between an interface system and two isolated surface slabs: (10) Here i is the surface energy of slab-i, 12 is the interface energy, i tot is the total energy of slabi, 12 tot is the total energy of a system (slab-1 + slab-2) with an interface, and A is the whole interface-area.

Results and discussions
Available analytical models for predicting the W sep were still generally limited to liquid-metal and oxide interfaces, and had to rely on some simple empirical correlations that described either the free energy of formation for oxide and liquid-metal, or the mixing enthalpies for respective oxide elements in metals [42].Unfortunately, many of these models were not applicable to those systems comprised of non-oxidized ceramics interface.Furthermore, these models may only provide qualitative information for trends in adhesion and their generality.
Even within the class of metal/oxide interfaces, it was still debatable because many of them only have those interfaces using the α-Al 2 O 3 (Alumina) as oxides [43].
In recognition of limitations of analytic models, a number of ab-initio studies on metal/ceramic adhesion using density functional theory (DFT) have appeared during past two decades.These studies captured atomic-scale features at the interfacial areas and can reveal invaluable information using electronic structures.For example, one DFT research focused on the interfacial adhesion at several Al-Ceramic couples: Al/α-Al 2 O 3 , Al/WC, Al/VC, Al/VN, Al/CrN and Al/TiN [44][45][46][47].Among which, large alloy slabs obeying a unit stacking sequence along x-and y-directions, respectively, but up to 15 atomic layers along z-direction were adopted in constructing calculation geometries, so as to accurately evaluate the surface energies and avoid the non-convergence problem via the slab thickness.Superscripts on polar surface systems, α-Al 2 O 3 (0001) O , WC(0001) W , WC(0001) C , etc., indicated surface terminations.In addition, each non-stoichiometric slab must be constructed to allow for the identical termination of its both surfaces.Regarding their relevant geometric parameters, see [47].In this research, although Eq. ( 10) defined the W sep in terms of both surface and interface energies, these quantities do not necessarily play an equal role, that is, since interface energy γ is relatively small in those metalceramic systems, the strength of interfacial bonding relevant to the W sep would depend to a large extent upon the reactivity of individual surfaces as reflected by their surface energies, see If considering atomic relaxations into the calculation, at an Al/Alumina interface, both magnitude and rank ordering of adhesion energies for different stacking sequences are found to change to the unrelaxed results determined by a "universal binding energy relation (UBER, in this framework, both charge density and binding energy of surface via inter-particle separation were empirically described as exponential decaying functions)" [48], see

Effect of surface geometry to the W sep
The W sep calculations are listed in Table 2 [47], where readers may find out obvious differences in the W sep values between polar and non-polar geometries; among which, three polar surfaces behave the largest W sep values from 4.08 J⋅m -2 for the Al/WC W up to 9.73 J⋅m -2 for the Al/α-Al 2 O 3 O ; for the remaining six non-polar interfaces, five have the W sep values of 2.10 J⋅m -2 or less.Furthermore, values for three nitride-ceramics are among the smallest range and look relatively insensitive to the choice of metallic component, with the W sep values falling within 1.45-1.73J⋅m -2 .Averagely, the W sep values in polar systems are about 4.60 J⋅m -2 larger than in nonpolar systems, consistent with their larger surface energies as shown in Table 1.It was also noted that, the predicted W sep value of 1.06 J⋅m -2 for the Al/α-Al 2 O 3 Al system was in favorable agreement with the experimental value of 1.13 J⋅m -2 in [49].

Effect of surface energies to the W sep
To further verify whether the W sep is mainly determined by surface energies or not, one Figure (see Figure 1 in [47]) plots the W sep versus the [σ Al +σ Ceramic ].A dashed-line with a unit slope is added to indicate data points lying in an idealized case: γ = 0.In this figure, the horizontal deviation of a given data point from the "γ = 0 -line" corresponds to the γ size.For three polar interfaces in the top region, data are bracketed with horizontal error bars, giving the range of possible σ i listed in Table 1.The W sep values in the whole range indicate a roughly linear trend with respect to the σ i , all but one of data points fall to the line-right.One system falling to the line-left with a negative γ value and involving an α-Al 2 O 3 O surface, shows the largest surface energy and the W sep value.This implies that the possibility of the fracture be within the Al-bulk rather than at the Al-Ceramic junction (interfacial area) due to a tensile stress.It also indicates that such an interface geometry is very firmly constructed.
With the exception of TiN, the rocksalt-structured carbides and nitrides (CrN, VN and VC) exhibit the largest dependence of the W sep on their surface energies: These three systems have the smallest interface energies, only 0.03-0.18J⋅m -2 , and their W sep values crowd in a nearly linear fashion (close to the γ = 0 -line) along the bottom left of the figure.Presumably, small γ values at these interfaces can be explained by similarities between bulk and interface bonding.However, since both of bulk components exhibit some extent to metallic bonding, for the Al/TiN system-a singularity with respect to other rocksalt ceramics, the γ value contributes larger, 0.62 J⋅m -2 , which could be due to the relatively large interfacial strain (misfit of slab geometry) of 5.3%.
Correspondingly, Table 2 also illustrates that those W sep values involving non-polar α-Al 2 O 3 and WC (1120) systems are not solely depicted by their surface energies (i.e., the γ effects should be more important).For two α-Al 2 O 3 interfaces, as with TiN, a possible explanation for this behavior could be the significant strain (4.9 and 5.3%).In addition, a main effect of ionic bonding in bulk α-Al 2 O 3 differs substantially from that in Al-bulk.On the contrary for the WC system, the strain is very low (0.4%), and thus differences of bonding types in respective bulks are less appeared (WC is metallic).So it implements an incoherent, misfit geometry characterized by irregular bonding across the interface.

Conclusion
Comprehensively, the DFT calculations indicated that the surface energy played a dominant role in determining the adhesion energy (the W sep ) for systems described above, interface energy was generally much smaller in comparison.In most cases, particularly for those polar systems involving cubic oxides, carbides and nitrides (adhesives) at interfaces, the W sep was described largely in terms of surface energy alone.Exceptions to this rule were mainly for systems with relatively large interfacial strain or incoherent interface bonding.

Ab-initio molecular dynamics study at an Al/Fe 2 O 3 interface
In this section, without any loss of generality, ab-initio molecular dynamics (AIMD) based upon density functional theory (DFT-GGA) is employed to discuss the adhesive transfer of materials at an interface combining a clean Al(111) and a hematite [α-Fe 2 O 3 (0001)] slabs together, to obtain invaluable information on the dynamical adhesion.As a result, the compression of two slabs may form an amorphous alumina (Al 2 O 3 )-texture, which welds two slabs together at the interfacial area through a thermite reaction.Subsequently, a tensile process results in the removal of amorphous alumina-texture from the melting Al-slab.The overall reaction is highly exothermic to increase the temperature of the system.

Modeling geometry and simulation procedures
An AIMD study focusing on a clean Al(111) slab here involves a super-cell geometry with a periodic 4-layer Al(111)-slab (16 ions per layer) in the x-and the y-directions [50].This orthorhombic super-cell has three orientated vectors: 110] = 2 2 Al at the x-axis; 112] = 6 Al at the y-axis; and 111] = 40 Å at the z-axis along with a vacuum distance of 24.0 Å in the z-direction, and one bottom layer fixed along the z-direction.Here Al = 4.053 Å [51].Similarly, simulations for a widely studied α-Fe 2 O 3 (0001) slab with oxygen-termination (more reactive than iron-termination to the Al-surface) are conducted using a periodic sevenlayer rectangle slab with two oriented vectors: 5 Fe 2 O 3 at the x-axis and Fe and 54 O ions).Here Fe 2 O 3 = 5.031 Å [52].Comparing these two slabs, the Al-slab size is slightly larger than the Fe 2 O 3 -slab along x-and y-directions, respectively.So a super-cell fitting for the Al-slab is selected as the final frame to contain both of slabs in the simulation.In this super-cell, misfit percentages of each of two slab lengths along x-and y-directions are calculated at very small values: 1.9 and 1.9%, respectively (see Figure 2).
The simulation time step was selected as 0.001 ps (i.e., 1 ps = 10 −12 s).The whole system was held at constant energy of AIMD simulation (no temperature control) when two slabs interacted.At the beginning of each simulation, top three layers of the α-Fe 2 O 3 -slab were constrained to move down along the z-direction at a constant velocity, V Z = -4.00Å/ps, corresponding to a displacement d Z = -0.004Å per time step.Since this speed, 400 m/s, was typical of real processing conditions, such shock loading may exceed the yield strength of the Al-slab, and the deformation of top layers in the Al-slab can be reasonably observed during contact simulation (Figure 2).Before carrying out the compressive simulation, a favorable interface geometry between two slabs is optimized using the energy minimization in the DFT.The result indicates that the final spacing between two slabs would be about 1.88 Å to form Al-O bonds at the interfacial area, which is in very agreement with those observed in experiments [53] (see Figure 3).Table 3 also lists some relevant bond lengths [50].[a] Ref. [55].

Adhesives -Applications and Properties
[b] Ref. [56].Therefore in the following AIMD simulation, the initial spacing between two slabs was selected as 2.50 Å, which was far enough to avoid the significant forming of chemical bonds in the R-L region (the interfacial area) when the simulation started.As two slabs moved closer, chemical bonds would begin forming in the spacing, and the interfacial reactions may occur.In      d) as an optimal bond spacing between two slabs is reached; and then, the stress increases again (after point e) when two slabs continue compressing.This figure also represents a polynomial (dashed line) fit to the stress curve, which may clearly regress the real stress trend through its oscillation.Please note that the stress trend just depicts the starting stage of a whole P-h curve as stated in [54].However, since two slabs here are too thin and too small in dimension, and at too short time steps throughout the whole AIMD simulation, it is hardly to rough out the whole P-h curve fulfilled by classical molecular dynamics on very large size models at very long time steps.Even if, this stress curve may still provide some invaluable information at the initial stage of nanoindentation process.

Tensile process
After the compression, an effect of a tensile strain acting on the compressed α-      Comprehensively, the AIMD simulations reveal that the tensile stress at the Al(111)/α-Fe 2 O 3 (0001) interface may remove the Al 2 O 3 fragments from the Al-slab surface after a steel roller (covered with an iron oxide) contacts the slab under an external strain load, which could indicate a typical process of the adhesive transfer of materials at interfaces.

Conclusion
The AIMD simulations of adhesive transfer of materials at an interfacial area without lubricants were performed for compressive and tensile processes between the Al(111) and the α- (3) The initial stage of the P-h curve can be drawn throughout the AIMD simulation to provide the invaluable information, which was hardly observed throughout classical molecular dynamics simulation.(4) During the AIMD simulation, the whole reaction process at the interfacial area was 2Al + 3Fe 2 O 3 → 6FeO + Al 2 O 3 , that is, due to the very short simulation time, such a thermite reaction was incomplete.

Summary
In this chapter, atomistic simulations of adhesive behaviors at several Al/Ceramic/Iron oxide interfaces are comprehensively discussed using density functional theory and ab-initio molecular dynamics, to illuminate a reasonable way for people's further understanding to this fascinating field.Some main points in the discussion include, without any loss of generality: 1. Main phenomenon in the Al/Ceramic adhesion indicated: In most cases, surface energy would play a dominant role in determining adhesion energy (W sep ) as interface energy was generally much smaller in comparison if ignoring effects of large stain and incoherent bonding at the interfacial area between two contact slabs, particularly for those polar systems involving cubic oxides, carbides and nitrides to form very solid thin films (adhesives) at interfaces.Apparently, the more the W sep depended upon the surface energy, the more steady and tenacious the adhesives (thin films) bind the interfaces.Therefore, knowledge of surface energies alone could serve as a reasonable starting point for estimating interfacial strength.

Figure 1 (
a). Regarding the optimized interfacial stacking sequences, through analyzing electronic structures (charge density contours), it is found that, at the polar Al/α-Al 2 O 3 and the polar Al/ WC interfacial areas (more dangling bonds on isolated surfaces), metal atoms prefer to the sites which continue the natural stacking sequence in ceramic bulk across the interfacial area into the Al-slab host, see the right picture in Figure 1(b); while for non-polar Al/α-Al 2 O 3 Al , Al/ VC, Al/VN, Al/CrN, Al/TiN and Al/WC systems, Al atoms prefer to the sits above metalloid atoms (O, C and N), see the left picture in Figure 1(b).

Figure 1 .
Figure 1.(a) (Left) An optimized geometry for a fcc-Al/α-Al 2 O 3 interface determined by the UBER calculations; (right) the final relaxed structure: Small spheres are Al atoms, and large spheres are O atoms.The view direction is along [1210]. (b) Side-view of electronic structures at an interface: Two slices through the electron localization function (ELF) for an fcc-Al/α-Al 2 O 3 interface taken along (1010) and (1120) planes, indicating four of hexagonal close- packed oxygen layers in oxide (bottom) and all five atomic layers from an Al slab (top).

15 2 Fe 2 O 3 at
the yaxis, and a height of 8.40 Å vacuum space with free boundary along the z-direction (total of 36

Figure 4 (
a), since top three layers (outside the R-L region) in the α-Fe 2 O 3 (0001)-slab may move down at a constant velocity (V Z = − 4.0 Å/ps) as a rigid block, its final downward displacement is chosen as 0.62 Å, so as to keep the final spacing between two slabs around an optimized value: 1.88 Å.

Figure 4
Figure 4 shows the compressing sequence of the Al(111) and the α-Fe 2 O 3 (0001) slabs.A fixed velocity, V Z = − 4.00Å/ps, is applied onto top three layers in the α-Fe 2 O 3 (0001)-slab along the z-direction during the simulation [see Figure 4(a)].

Figure 5 .
Figure 5.The stress perpendicular to the interfacial area as two slabs compress together.

Figure 5
Figure5indicates the distribution of the stress normal to the interfacial area as two slabs compress together.Initially, the stress becomes negative (attractive between two slabs) as Al-O bonds begin forming in the R-L region (points a-c); when two slabs are moving closer, the stress tends to the positive values (repulsive between two slabs) because most of Al and O atoms at the interfacial area are forced to come closely, leading to a nearly formation of amorphous alumina-texture with thermal (exothermic) expansion; after two slabs totally react, the stress falls down to a near zero point (point e) from the highest positive one (point d) as an optimal bond spacing between two slabs is reached; and then, the stress increases again (after point e) when two slabs continue compressing.This figure also represents a polynomial (dashed line) fit to the stress curve, which may clearly regress the real stress trend through its oscillation.Please note that the stress trend just depicts the starting stage of a whole P-h curve as stated in[54].However, since two slabs here are too thin and too small in dimension, and at too short time steps throughout the whole AIMD simulation, it is hardly to rough out the whole P-h curve fulfilled by classical molecular dynamics on very large size models at very long time steps.Even if, this stress curve may still provide some invaluable information at the initial stage of nanoindentation process.

Figure 6
Figure 6 shows the number change of Al-O and Fe-O bonds versus normal distance at the interfacial area in Figure 4(a) during the compressive process.In this figure, Al-O bonds in the R-L region increase substantially while Fe-O bonds in this region decrease gradually.The reason is that, when two slabs compress together, the dangling bonds on Al and O atoms at the interfacial area would bond together.As the compression continues, more Al-O bonds are forming as some Fe-O bonds are breaking.Nevertheless, only a very porous interface is formed during two slabs compression, which is far from the ideally compressed status.Therefore, the thermite reaction at the interfacial area just partially occurs, only a few Fe-O bonds are breaking in this region.

Figure 6 .
Figure 6.The change of Al-O and Fe-O bonds in the R-L region during slab compression.
Fe 2 O 3 (0001)/ Al(111) slabs is immediately conducted: Top three layers of the α-Fe 2 O 3 (0001)-slab in Figure 4(e) is moving up as a rigid block at a constant velocity, V Z = + 4.00 Å/ps, reversing the downward direction of previous compression.

Figure 7
Figure7shows the stretching sequence of this process: (i) when top three layers of the α-Fe 2 O 3 (0001)-slab move up about 1.30 Å from the starting point in Figure4(e), the real tensile process at the interfacial area begins (part a).(ii) Then next, as two slabs separate, further deformation at the interfacial area and thermite reaction will continue, more atomic mixing in this area oxidizes the Al-surface (part b).(iii) During transitional stages (parts c and d), the whole interfacial area begins breaking apart, mainly due to the breaking of some more Fe-O and Al-Al bonds.Along with a small increase of Al-O bonds, amorphous alumina-texture is continuously formed, and an adhesive metal transfer begins.(iv) Finally, surface on the Alslab is totally broken: an adhesive metal transfer has completely occurred, amorphous alumina-texture is removed from the remnant surface of the Al-slab, and bonds onto the

Figure 8 .
Figure 8.The number change of Al-O, Fe-O and Al-Al bonds in the new R-L region during the slab tension.

Figure 8
Figure8shows the number change of Al-O, Fe-O and Al-Al bonds in a new larger R-L region as shown in Figure7(a) during the tensile process.In this figure, at initial stage, a few more Al-O bonds are forming and a few more Fe-O bonds are breaking, which represents the continuous thermite reaction in the R-L region as the Fe 2 O 3 -slab is being pulled upward.Also, there is a breaking of approximately 10-12 Al-Al bonds in top two layers of the Al-slab.When the Fe 2 O 3 -slab is pulled upward about 8.0 Å, there is a slight increase in Al-O bonds but a slight decrease in Fe-O bonds.In this case, the interfacial area is totally separated.By the end of the tensile process, the Fe-O coordination is approximately 1:1, while the Al-O one is approximately 2:3.Thus, this reaction could be approximately 2Al + 3Fe 2 O 3 → 6FeO + Al 2 O 3 , instead of a more favorable thermite reaction: 2Al + 3FeO → 2Al 2 O 3 + 3Fe, which is mainly due to very limited reaction time.Local temperature in the residual Al-slab is about 810 K, while that in the removed slab (Fe 2 O 3 /FeO+Al 2 O 3 ) is about 3030 K, implying that such a soft slab may bond onto a hard-roller surface as a protective thin film resistant to the elevated temperature.

Fe 2 O 3 (
0001) slabs.The conclusions can be drawn: (1) Contact of the Fe 2 O 3 -slab with the Al-slab may result in a strong exothermic reaction to form an amorphous alumina-texture, which welded two slabs together.(2) Separation of two slabs may have the local temperature at the interfacial area increased significantly, leading to the melting of the Al-surface and easy removal of amorphous Al 2 O 3 -texture from the Al-slab surface.

2 .
Common phenomena in the adhesive transfer of materials (phase transformation of indented materials) included: (a) Contact of the Fe 2 O 3 -slab with the Al-slab may result in a strong exothermic reaction to form an amorphous alumina-texture, which welded two slabs together.(b) Separation of two slabs led to the melting of the Al-surface and easy removal of amorphous Al 2 O 3 -texture from the Al-slab surface.This texture may bond onto a hard-roller surface as a protective thin film (soft adhesives) resistant to the elevated temperature.

Table 1 .
Calculated surface energies σ i and slab-i dimensions (no. of layers) used for the W sep calculations.
In this figure, at initial stage, a few more Al-O bonds are forming and a few more Fe-O bonds are breaking, which represents the continuous thermite reaction in the R-L region as the Fe 2 O 3 -slab is being pulled upward.Also, there is a breaking of approximately 10-12 Al-Al bonds in top two layers of the Al-slab.When the Fe 2 O 3 -slab is pulled upward about 8.0 Å, there is a slight increase in Al-O bonds but a slight decrease in Fe-O bonds.In this case, the interfacial area is totally separated.By the end of the tensile process, the Fe-O coordination is approximately 1:1, while the Al-O one is approximately 2:3.Thus, this reaction could be approximately 2Al + 3Fe 2 O 3 → 6FeO + Al 2 O 3 , instead of a more favorable thermite reaction: 2Al + 3FeO → 2Al 2 O 3 + 3Fe, which is mainly due to very limited reaction time.Local temperature in the residual Al-slab is about 810 K, while that in the removed slab (Fe 2 O 3 /FeO+Al 2 O 3 ) is about 3030 K, implying that such a soft slab may bond onto a hard-roller surface as a protective thin film resistant to the elevated temperature.