A. Physical models and computational methods

To be able to study, in oxygen environment, the evolution of metallic surface steps created by the emergence of dislocations and the formation of micro-notch during cyclic loading, we use sophisticated phenomenological models such as reactive force field (ReaxFF) and/or charge-optimized many-body potential (COMB) as implemented within the large-scale atomic/molecular massively parallel simulator (LAMMPS) code. ^{Reference Plimpton26,Reference Plimpton, Thompson, Shan, Moore and Kohlmeyer27} These models were developed to simulate interactions between different types of elements such as nickel and oxygen or copper and oxygen. Reliable simulation of interaction between metallic atoms and oxygen molecules or oxygen anions implies that each atom can dynamically and autonomously determine its charge according to its local environment. This can be achieved by a self-consistent charge equilibration scheme based on the principle of electronegativity equalization (EE) thanks to either the EE method (EEM) developed by Mortier et al. in 1986, ^{Reference Mortier, Ghosh and Shankar28} or the electronegativity equilibration (QEq) method developed by Rappe and Goddard III in 1991. ^{Reference Rappe and Goddard29} ReaxFF uses EEM, ^{Reference Senftle, Hong, Islam, Kylasa, Zheng, Shin, Junkermeier, Engel-Herbert, Janik, Aktulga, Verstraelen, Grama and van Duin30} whereas COMB uses QEq coupled with Tersoff-like potentials. ^{Reference Shan, Devine, Kemper, Sinnott and Phillpot31,Reference Liang, Shan, Cheng, Devine, Noordhoek, Li, Lu, Phillpot and Sinnott32} Ab initio calculations are unfortunately clearly ruled out for the study at stake because of time and memory-space limitations of currently available computers.

ReaxFF sets developed for nickel and oxygen ^{Reference Assowe Dabar33} and for copper and oxygen ^{Reference Psofogiannakis, McCleerey, Jaramillo and van Duin34,Reference Van Duin, Bryantsev, Diallo, Goddard, Rahaman, Doren, Raymand and Hermansson35} will be used. COMB3 will also be used for copper and oxygen. ^{Reference Liang, Shan, Cheng, Devine, Noordhoek, Li, Lu, Phillpot and Sinnott32} On top of the tests of these potentials which were already carried out, ^{Reference Liang, Shan, Cheng, Devine, Noordhoek, Li, Lu, Phillpot and Sinnott32–Reference Van Duin, Bryantsev, Diallo, Goddard, Rahaman, Doren, Raymand and Hermansson35} we examine the stiffness elastic constants of both metals, Ni and Cu, and their cubic oxide counterparts. The main goal is to compare these stiffnesses from an intrinsic point of view for Ni and Cu. Yet, it is important to keep in mind that hereafter in Sec. III we will not simulate clean interfaces between metals and oxides but the fatigue of metals in oxygen environment. The study of the structures of the early stages of surface oxidation, starting with oxygen deposition and the formation of specific two-dimensional substructures, will also be carried out in Sec. III.

The cubic lattice parameters and elastic constants are gathered in Table I for Ni and NiO [nickel monoxide or nickel(II) oxide] as well as for Cu and Cu₂O (cuprous oxide). Cupric oxide, CuO, is monoclinic. ReaxFF current descriptions of Ni and Cu (compatible with O) are restricted to a simple pair potential description which forces the Cauchy relation c ₁₂ = c ₄₄. ReaxFF has been primarily developed to study superficial chemical reactions and not bulk properties. This ought to be improved in the future. COMB3 is more successful to describe copper. True NiO is magnetic and rhombohedral below its 523-K Néel temperature. Yet, it is almost cubic, with a room temperature rhombohedral angle α_r ∼ 60.4° very close to the corresponding cubic value of 60°, and we can consider it as cubic. The elastic constants c ₁₂ and c ₄₄ are numerically difficult to estimate with ReaxFF, and a negative value was also obtained by Assowe Dabar ^{Reference Assowe Dabar33} for c ₁₂ which fortunately does not correspond to a deformation mode by itself. It should be noted that experimental values themselves show very large deviations (and the experimental values from Ref. Reference Wiame, Maurice and Marcus46 actually corresponds to weighted values from four sets of values showing deviations). Ab initio calculations, à la Hartree–Fock or within the local spin density approximation plus a Hubbard term also exhibit large deviations. ^{Reference Towler, Allan, Harrison, Saunders, Macrodt and Aprà37,Reference Dudarev, Botton, Savrasov, Szotek, Temmerman and Sutton38} Comparing the stiffness elastic constants of Ni and NiO, it is difficult to say which one, Ni or NiO, is stiffer. The two available sets of experimental stiffness constants for cuprous oxide are more consistent than for nickel oxide but are not reproduced too well by COMB3 (and are difficult to estimate with ReaxFF). Yet, for Cu and Cu₂O, one can say that Cu is stiffer than Cu₂O. Although we will not be dealing with either bulk NiO or bulk Cu₂O per se in the following, we think it is interesting to note that none of these two oxides are stiffer than their pure metal counterparts. The melting temperature T _m of Ni predicted by ReaxFF potential is about 1950 K, not too far from the experimental T _m which is 1728 K. The T _m values predicted by ReaxFF and COMB for copper are about 1450 and 1300 K respectively, reasonable in comparison to the experimental value 1358 K. As these potential-dependent T _m values are not too far from the experimental T _m, we shall indicate our simulated temperatures T numerically in Kelvin rather than in reduced units (T/T _m). Other tests of the efficiency of these potentials can be found in the articles where they have been developed. ^{Reference Liang, Shan, Cheng, Devine, Noordhoek, Li, Lu, Phillpot and Sinnott32–Reference Wang, Fisher and Manghnzmi36} No modeling being perfect, we intentionally used several potentials to ensure the validity of our final conclusions.

TABLE I.

Lattice parameter a (in Å) and stiffness elastic constants (in GPa) of Ni, NiO, Cu, and Cu₂O. The two cubic shear elastic constants are c ₄₄ and c′ with c′ = (c ₁₁ − c ₁₂)/2. The isotropic incompressibility is given by the bulk modulus B = (c ₁₁ + 2c ₁₂)/3. Our values are calculated at 0 K. Experimental values are given as extrapolated at 0 K (Exptal 0 K) and at 300 K, viz. room temperature (Exptal RT). The Ni and NiO ReaxFF values calculated in this work are similar to the ones calculated in Basinski and Basinski. ^{Reference Basinski and Basinski23} See main text for more comments.

^a Ref. Reference Simmons and Wang39.

^b Ref. Reference Du Plessis, van Tonder and Alberts40.

^c Ref. Reference Uchida and Saito41.

^d Ref. Reference Wiame, Maurice and Marcus46.

^e Ref. Reference Hallberg and Hanson42.

^f Ref. Reference Beg and Shapiro43.

The equilibrium distance between the two oxygen atoms of an O₂ molecule is about 1.2 Å, correctly reproduced by both ReaxFF and COMB. Neither ReaxFF nor COMB includes the magnetic properties of O₂ molecules which can indeed be neglected except at very low temperatures. ^{Reference Etters and Hardouin Duparc44,Reference Hardouin Duparc and Etters45} ReaxFF and COMB are phenomenological descriptive models, not quantum mechanical models. They allow us to consider larger systems on longer MD simulation times than what first principles calculation can allow. Simulations will be carried out at constant temperatures via the Nosé–Hoover thermostat implemented in LAMMPS, with a Verlet time step of 0.5 fs, small enough to ensure the stability of the algorithm. Two temperatures will be considered, 300 and 600 K.

B. Simulation systems and loading methods

In this work, we will first investigate the first oxidation stages of nickel and copper surfaces, specifically their {111} and {100} surfaces. Our results will be compared with experimental observations carried out on the same surfaces. This will allow us to both check the reliability of the ReaxFF and COMB models and to gain insights into the oxidation mechanisms at the atomistic level. We will then repeat in oxygen environment our previous study of surface slip irreversibility and nano notch initiation carried out in inert environment in Ref. Reference Fan, Hardouin Duparc and Sauzay5 in order to check the influence of an oxygen atmosphere. Three simulations systems will be used as described now.

The first simulation system [Fig. 1(a)] is built for simulation of oxidation of a {111} surface, i.e., the most compact surface for fcc metals. The simulation box configuration is presented in Fig. 1. The x, y, and z axes correspond to the $\left[ {1\bar{1}0} \right]$ , [111], and $\left[ {11\bar{2}} \right]$ crystallographic directions, respectively. The lattice periodicities in these directions being ${a_x} = a{\rm{/}}\sqrt 2$ , ${a_y} = a\sqrt 3$ , and ${a_z} = a\sqrt {3/2}$ , the dimensions of the metal crystallite are 14a _x × 4a _y × 14a _z , with about 5000 atoms of nickel or copper (14 × 4 × 14 × 6 atoms). Although bulk periodic boundary conditions are applied to directions x and z, a vacuum space about 120 Å (12 nm) large is introduced in the [111] direction, to simulate two (111) surfaces and to introduce oxygen molecules. About 2000 molecules are initially placed in that vacuum space with approximately one molecule per 150 ų. It corresponds to a very large gas density and will induce a high gas pressure. It is necessary to make the first steps of oxidation attainable at MD scales. The closest oxygen molecules to the (111) surfaces are initially about 20 Å farther from them, so that these oxygen molecules have time (a few picoseconds) to relax before they get in contact with the metal surfaces. The diagonal stress components along the metal bulk periodic directions are controlled within the LAMMPS code, set equal to about zero in average. The stress component along the third direction corresponds to the oxygen pressure. It is first estimated to be about 6 × 10⁷ Pa at 600 K in the perfect gas approximation, and then kept constant in average. We recall that imposing stresses means that the box sizes are adaptative.

FIG. 1.

x-y plane views of the three simulation boxes used in this work. Red atoms are metal atoms and blue atoms are oxygen atoms. The right third simulation box is presented by a sketch in order to show clearly the slip system. The two left boxes correspond to snapshots at 0.5 ps after the start of the MD simulation.

The second simulation system [Fig. 1(b)] is aimed at analyzing the oxidation of the open {100} surface. The x, y, and z axes are oriented according to the [100], [010], and [001] directions, respectively. The dimensions of the metal substrate are 14a × 6a × 14a with a total of 4704 atoms. As in the previous case, the boundary conditions are applied with a vacuum space of 125 Å in the [010] direction initially filled with about 1400 oxygen molecules. We use the same temperature and pressure conditions.

The third simulation system [Fig. 1(c)] is used to analyze surface slip irreversibility in oxygen environment and to compare with simulations previously done in inert environment. ^{Reference Fan, Hardouin Duparc and Sauzay5} We have a periodically repeated metal bar which will be surrounded by oxygen molecules. The three axes are x: $\left[ {\bar{1}31} \right]$ , y: [714], and z: $\left[ {11\bar{2}} \right]$ . The (111) slip plane is perpendicular to the x–y plane and the $\left[ {1\bar{1}0} \right]$ slip direction belongs to the x–y plane. Perfect b = $a\left[ {1\bar{1}0} \right]/2$ edge dislocations are introduced by inserting two half-planes in the bulk of the simulation box; the other atoms being displaced using the simple isotropic elastic approximation. The system is then relaxed with MD with the chosen Ni or Cu interatomic potential and perfect dislocations dissociate into two Shockley partials. The dislocations can glide to the $\left( {\bar{1}31} \right)$ surface under small shear stress. The typical dimensions of the simulated metal box in this work are 44, 266, and 24 atomic planes along x, y, and z, respectively, roughly corresponding to 4.9, 6, and 1.75 nm for copper with ReaxFF. There are nearly 4300 metal atoms. Casual simulations with larger systems and more atoms showed that it does not influence our results. Bulk periodic boundary condition is used along the z-direction with a length period adapted to take into account the thermal expansion. A 60-nm large vacuum space is introduced in the x direction, initially filled with almost 500 oxygen molecules. Simulations have been carried out at both 300 and 600 K. The oxygen pressure turns out to be nearly 5 × 10⁷ Pa at 600 K. To study slip irreversibility with oxidized $\left( {1\bar{3}1} \right)$ surfaces, we first let oxidation occur during 2 ns at 600 K to get a thick-enough oxide layer. The system is then cooled down to room temperature (RT) and dislocations are inserted in the bulk and glided to the oxidized surface located on the right.

The upper and lower horizontal boundaries are monitored to initially apply a uniform tensile or compression displacement in the y direction. These upper and lower monitored walls are typically 3-Å thick [which roughly corresponds to thirteen (714) “atomic planes,” although these planes are not well distinguishable atomic planes as the (111) atomic planes are]. Atoms within these y-monitored layers are free to move along the x and z directions. In the case of cyclic loading, repetitive continuous sinusoidal elongation–compression displacement is applied to the upper and lower boundaries along the y direction. The strain rates are in the order of 10⁷–10⁸ s⁻¹, typical in similar MD simulations.

Just as strain rates are taken much larger in MD simulations than in actual experimental fatigue tests to counterbalance the smallness of the simulation times currently available for systems containing thousands of atoms, oxygen pressure in our MD boxes is much larger than the actual oxygen pressure in atmospheric conditions to be able to get some oxidation stages during accessible MD simulation times. A corresponding measure of the exposure of metal to oxygen is the product (oxygen pressure) × (time of exposition). The standard unit of exposure traditionally used by experimentalists is the Langmuir L unit, with 1 L = 1.33322 10⁻⁴ Pa s (because the original definition is 1 L = 10⁻⁶ Torr s). In our simulations at 300 K, with a pressure of 3 × 10⁷ Pa during 2 ns, the exposure is 450 L. This oxygen exposure is comparable to the one achieved in several experimental studies of first stages oxidations: for instance, ^{Reference Wiame, Maurice and Marcus46} has oxygen exposure of 900 L at RT. For a very high cycle fatigue (VHCF) test of frequency of 10,000 Hz at natural atmospheric conditions (hence a partial oxygen pressure of about 0.2 × 10⁵ Pa), the exposure is typically 15,000 L during one cycle. This exposure is about thirty times larger than in our MD simulation at 300 K, and fifteen times at 600 K. However, our simulations showed that the chemisorption of oxygen occurs very rapidly at the beginning of the simulation and becomes much slower after typically 1 ns at 600 K (Fig. 2). Thus, proportionally, only few more O₂ molecules are chemisorbed at 15,000 L than at a hundred of langmuirs so that our MD simulations can be thought to be meaningful with respect to, at least, the VHCF tests.

FIG. 2.

Number of chemisorbed oxygen atoms at the $\left( {14a\sqrt 2 } \right)$ × $\left( {14a\sqrt {\left( {3/2} \right)} } \right)$ nickel (111) surface as a function of time, at 600 K and 6 × 10⁷ Pa. Per conventional definition, the (hkl)-monolayer coverage is equal to the (hkl) atomic planar density of the substrate. For a (111) surface, it is $4{\rm{/}}\sqrt 3$ atoms per a ², hence 392 atoms for our surface. Of course, oxygen atoms do not permanently stay on the surface and penetrate the nickel substrate and oxidize it, as described in the main text and as shown in Fig. 3.

A. Nickel

1. First stages of nickel surface oxidation

Oxidation is first studied on nickel (111) and (010) surfaces to compare it with what is experimentally known on these surfaces. Simulations are carried out at 300 and 600 K. The same interaction mechanism, to be described below, is found between oxygen and nickel at these two temperatures.

According to experimental investigations, ^{Reference Fehlner and Mott47–Reference Holloway50} three reaction regimes can be identified during nickel oxidation: a very rapid chemisorption, nucleation of NiO islands at the surface and their coalescence, and a stable growth of the oxide layer. The first regime and the beginning of the second regime are observed in our simulations, providing us important atomistic details.

a. First oxidation stages of the {111} nickel surface

Oxidation is first studied on nickel (111) surface. We found that once an oxygen molecule gets on the nickel surface, it dissociates into two oxygen atoms very rapidly, typically within 5 femtoseconds. No stable state of chemisorbed oxygen molecule is observed, in agreement with the experimental studies mentioned above. The number of dissociated chemisorbed oxygen atoms increases very fast during the first 10 ps, and then the adsorption rate decreases significantly (Fig. 2). The adsorption rate depends on the number of particles striking the surface per second and on the so-called “sticking coefficient,” which is the probability that an impinging particle actually sticks to the substrate, at a place where there is no previously adsorbed oxygen. With the augmentation of the number of absorbed oxygen molecules on the surface, new incoming oxygen molecules tend to bounce back and the sticking coefficient drops down.

Adsorption on the (111) surface becomes more difficult once the oxygen coverage reaches nearly 0.35 (111)-monolayer (ML). This suggests a correspondence with the formation of the 1/3 ML (111) $\left( {\sqrt 3 \times \sqrt 3 } \right)$ R30° superstructure (Fig. 3), following Wood’s classical notations. ^{Reference Wood51} The superstructure becomes more and more regular during the simulation. Lower oxygen coverage superstructure such as the (111) (2 × 2) ^{Reference Davisson and Germer52} has not been observed at the beginning of the oxidation dynamics, most certainly because our pressure and temperature conditions do not allow this superstructure to develop.

FIG. 3.

Oxidation on the Ni(111) surface at 600 K and 6 × 10⁷ Pa. Oxygen atoms are shown in blue and nickel atoms are in red [and yellow in (a₂), see below]. Two zones with $\left( {\sqrt 3 \times \sqrt 3 } \right)R30^\circ$ structures can be identified in (a₁), a snapshot taken at 255 ps. Three (111) nickel atomic layers are shown in (a₂), corresponding to the (111) ABC periodic stacking. (a₂) is an enlarged view of an instantaneous Ni(111) $\left( {\sqrt 3 \times \sqrt 3 } \right)R30^\circ$ oxygen superstructure with oxygen atoms being more or less in hollow sites. The light red and yellow atoms are the topmost nickel atoms. The nickel atoms that are on the layer below are in dark red and the third layer of atoms are in purple. Some top-most nickel atoms are colored in yellow to better show the superstructure. (b) A side view showing oxygen penetration, nickel diffusion, and the formation of an “external layer” on the top of the topmost layer of the original nickel substrate after 2 ns. The oxygen atoms are between the nickel (111) layers. The structure of the corresponding external layers is given in (c), showing two quasi NiO bidimensional (2d) germs.

Oxygen atoms diffuse into the bulk. Some oxygen atoms penetrate until the fourth atom layer after 2 ns at 600 K and 6 × 10⁷ Pa. The metallic slabs being twelve atomic layers thick, they are large enough for the simulations carried out in this study. Oxygen atoms are located in octahedral interstitial sites between the (111) atomic planes (the interatomic plane distance is initially $a{\rm{/}}\sqrt 3$ ) [Fig. 3(b)]. By locally filling the octahedral sites between the nickel atoms with oxygen atoms, the tendency is to progressively build the stable NiO compound (space group $Fm\bar{3}m$ , #225, Pearson symbol cF8, fcc with a specific two-atom motif per node, the only unambiguous symbolic description being the Strukturbericht type symbol B1, or two interlocked fcc elementary structures, each one with atoms being in the octahedral sites of the other). Yet the lateral distances remain constrained by the nickel bulk substrate. The nickel atoms which remain in the substrate do not seem to be highly disturbed, from a positional point of view, by the additional presence of oxygen atoms in the interstitial sites. They keep their fcc relationship. Whereas some oxygen atoms have diffused inside the nickel substrate, some nickel atoms have conversely diffused outward, inducing an external growth, with an external layer emphasized in Fig. 3(b). The nickel atoms on the new, external layer are intermixed with the oxygen atoms and form two-dimensional islands of NiO structure as can be seen in Fig. 3(c) which is a snapshot obtained after a simulation time of about 2 ns. Within these islands, the nearest nickel–nickel distance is nearly 3 Å, significantly larger than the nearest 2.55 Å distance in the fcc substrate and typical of the nearest distance between the nickel atoms in a perfect cubic (B1 or NaCl-like as mentioned above) NiO structure at 600 K (with a cubic lattice parameter nearly equal to 4.2 Å).

b. First oxidation stages of the {001} nickel surface

As for the {111} surface, there is a rapid dissociative chemisorption of oxygen molecules on the (010) surface, within the first picoseconds, followed by a significant decrease of the adsorption rate. A (010) c(2 × 2) superstructure is observed (Fig. 4), which corresponds to a (010) surface oxygen coverage of 50%. Oxygen adsorption is more difficult afterward. Similar to the (111) surface case, lower oxygen coverage superstructure (010) p(2 × 2) has not been observed at the beginning of the oxidation dynamics, due to applied pressure and temperature conditions.

FIG. 4.

Oxidation on the 14a × 14a Ni(010) surface at 600 K and 6 × 10⁷ Pa after 2 ns. (a) A side view showing oxygen penetration and the outward diffusion of a few nickel atoms. (b) The topmost nickel layer with oxygen atoms on top of it. The green rectangle embodies a {100} c(2 × 2) superstructure.

Oxygen atoms also diffuse into the bulk and penetrate until the fourth atom layer. Oxygen atoms are also located in octahedral sites and thus within the (010) planes, see Fig. 4(a) where the oxygen atoms seem to be less numerous only because they are hidden by the nickel atoms. Only a few atoms diffuse outward, even after a simulation time of about 2 ns, and no 2d NiO island is observed, as contrary to the (111) case. This difference can be ascribed to the fact that in the (010) case, oxygen atoms can relatively easily make their way through the substrate by diffusing between the nickel {111} atomic planes which make a 30° angle with the (010) surface. In the (111) surface case, the oxygen atoms have to considerably disturb the positions of some nickel atoms to diffuse through the (111) planes, allowing or even inducing the outward diffusions of some of these nickel atoms.

Simulations occasionally performed with larger systems showed that the dimension of our simulation box is large enough to investigate early stages of oxidation. Increasing the thickness does not bring new information. The formation of truly crystalline oxides is beyond the scope of the current simulations.

Our studies about oxidation on the {111} and {100} surfaces are in good agreement with existing experimental studies and bring more atomistic details. This preliminary study also indicates that the potential we use is likely to be good enough to be used for simulations of oxidation on the {131} surface, for which there are no experimental results.

c. First oxidation stages of the {113} nickel surface

Oxidation mechanisms on nickel {113} surface are similar to those on nickel {111} and {001} surfaces. There is a dissociative chemisorption of oxygen molecules, some oxygen atoms diffuse inward while some nickel atoms diffuse outward, generating an external layer with the formation of 2d NiO embryos again [Figs. 5(a) and 5(c)]. This time we chose to plot the nickel and oxygen atom densities along the direction perpendicular to the surface [Fig. 5(b)]. The small “peaks” at the right end correspond to the atoms in the external layer, and they correspond to a 1:1 stoichiometry typical again of the (100) plane of nickel monoxide NiO or of 2d islands of that structure.

FIG. 5.

Oxidation on the Ni $\left( {\bar{1}31} \right)$ surface at 600 K and 6 × 10⁷ Pa after 2 ns. (a) A side view showing oxygen penetration and nickel outward diffusion. (b) Distribution of nickel and oxygen atoms along the x $\left[ {\bar{1}31} \right]$ direction. To draw this atom distribution, the simulation box is cut by slices of 0.2 Å in the x-direction. The number of atoms N is counted in every slice and presented by square symbols which are linked to “show” the atom layers: every “peak” corresponds to an atom layer. The {113} atomic planes are relatively sparse (low atomic planar density) with small distances between them, viz. $a{\rm{/}}\sqrt {11}$ . (c) The external layer with the formation of small NiO islands.

The positions of nickel atoms at the initial surface plane are only slightly disturbed by the chemisorption of oxygen atoms. The explanation is the same as in the {001} surface case because the {111} nickel atomic planes make a 31.5° angle with $\left( {\bar{1}31} \right)$ surface plane.

2. Surface slip irreversibility of nickel in oxygen environment

Surface relief evolutions are now analyzed depending on the number of applied cycles. MD cyclic loadings are applied by repetitive continuous sinusoidal elongation–compression displacement along the y-direction, as described in the method Sec. II.B.

Surface step reconstruction is first tested in oxygen environment at five different temperatures, viz. 10, 200, 300, 500, and 700 K. An initial configuration with a 2b large surface step is geometrically prepared by cutting atoms away [Fig. 6(a)]. It is observed that the surface step is reconstructed at all temperatures in oxygen environment, as it is in inert environment. ^{Reference Fan, Hardouin Duparc and Sauzay5} The reconstructed surface steps are equally totally irreversible. 300 and 700 K test results are shown in Fig. 6. Snapshots are taken at each temperature after 30 ps. The oxidation and surface reconstruction take place at the same time, and the surface step reconstruction is even slightly accelerated, probably because oxygen atoms disturb nickel atom positions, which facilitates the step reconstruction. From these snapshots, one also sees that oxygen penetrates deeper with the augmentation of temperature and of oxygen gas pressure.

FIG. 6.

(a–c) Surface step reconstruction at 300 and 700 K after 30 ps in nickel with oxygen. It is similar to step reconstructions observed in inert environment. ^{Reference Fan, Hardouin Duparc and Sauzay5}

These simulations have been carried out with geometrically prepared surface steps to ease the visualization of step reconstruction. In a real fatigue condition, steps are created by dislocations gliding toward an already oxidized surface. In the following, we first let oxidation occur during 2 ns at 600 K on the $\left( {\bar{1}31} \right)$ surfaces. The system is then cooled down to RT, and dislocations are inserted in the bulk and glided to the oxidized surface located on the right.

Dislocations are introduced in the bulk (see method Sec. II.B) near the free surface as shown in Fig. 7. They can glide to the free surface under an applied small stress (∼1 MPa) [Figs. 7(a) and (b)]. A surface step is created and gets reconstructed afterward, similar to what was observed in inert environment. ^{Reference Fan, Hardouin Duparc and Sauzay5} Spontaneously or through thermal agitation, the tip atoms of the surface step, whether nickel or oxygen, move closer to the surface, optimize their cohesive energies by getting more neighbors, and thus minimize the initial surface step energy.

FIG. 7.

(a) and (b) Gliding of a dislocation under compression toward the oxidized surface, creation, and reconstruction of surface steps, at 600 K. The green lines represent dislocations in a purely schematic way. (c) The initiation of a surface micro-notch at 300 K in oxygen environment. Six opposite sign dislocations are inserted in the bulk in not direct neighboring (111) planes for the initiation of this micro-notch [see Fan et al. ^{Reference Fan, Hardouin Duparc and Sauzay5} The (111) planes of insertion are separated by four (111) atomic planes]. The oxide layer formed at 300 K is thinner than at 600 K.

Reconstructed steps are irreversible under reversed loading. A reconstructed step is reversible only if an opposite sign dislocation whose core is located in a direct neighboring plane arrives to the same surface. If there are one or several atomic planes between opposite sign dislocations arriving at the surface, the step is irreversible and a micro-notch forms under cyclic loading, as in inert environment [Fig. 7(c)]. The simulations for Figs. 7(a) and 7(b) have been carried out at 600 K so as to get relatively thick “oxide” layers to check their effects on step reconstruction. Since they have no significant effect, the simulation for Fig. 7(c) was carried out at room temperature.

B. Copper

1. First stages of copper surface oxidation

As mentioned in Sec. II.A, two atomic potentials exist for copper with oxygen, viz. COMB3 and ReaxFF. We will use both of them, at different temperatures. Results will be compared with experimental observations, ab initio simulations, and the previous results about nickel oxidation. The simulation systems are the same as for nickel oxidation simulations. We will mainly illustrate our results in this section on the $\left( {\bar{1}31} \right)$ copper surface because, as for nickel, the oxidation mechanisms for the three {111}, {001}, and {113} surfaces are essentially similar except for the differences which are due to the different {hhl} crystallographies.

Copper surface oxidation is different from that of nickel. In contrast to nickel oxidation, where oxygen molecule dissociation is almost immediate even at a relatively low temperature, viz. RT, we mainly observed chemisorbed states of oxygen molecules on copper with both potentials (Fig. 8). This is in agreement with experimental ^{Reference Wiame, Maurice and Marcus46,Reference Spitzer and Lüth53,Reference Sueyoshi, Sasaki and Iwasawa54} and ab initio studies. ^{Reference Lou, Nordlander and Hellsing55} In Devine et al., ^{Reference Devine, Shan, Cheng, McGaughey, Lee, Phillpot and Sinnott56} simulations with COMB3 potential showed that on a (001) surface, oxygen molecules dissociate only if they are on “top sites” above a prominent copper atom of the surface, and that if an O₂ is adsorbed on a “bridge site” or in a “hollow site” between surface atoms, it remains as an O₂ molecule. We found that with both COMB3 and ReaxFF, oxygen molecules tend to avoid going on top of prominent copper atoms and get chemisorbed at bridge or hollow sites. Chemisorbed oxygen molecules may finally dissociate, but it presumably takes a much longer time {until 2 ns of our MD simulations, COMB3 does not induce any oxygen dissociation when ReaxFF predicts some dissociations [Fig. 8(d)]}.

FIG. 8.

Chemisorption of oxygen molecules at a $\left( {\bar{1}31} \right)$ copper surface at 300 K after 50 ps, with COMB3 (a, b, c) and with ReaxFF (d, e, f) [two dissociated oxygen atoms are visible in (d)]. (a, d) are $\left[11\bar2} \right]$ views and (b, e) are $\left[\bar131} \right]$ views. (c, f) are equivalent $\left[\bar131} \right]$ views of only the top atoms for clarity, taken at a little earlier time than b and e. Simulations at 600 K give similar results, with more oxygen dissociations at 600 K with ReaxFF. (c) and (f) Oxygen molecules adsorption sites on the $\left( {\bar{1}31} \right)$ copper surface, as differently predicted by COMB3 and ReaxFF; oxygen molecules are mainly adsorbed on bridge sites with COMB3, but on hollow sites with ReaxFF.

We observed an important difference between the predictions made by COMB3 and ReaxFF. With COMB3, several oxygen molecules get adsorbed in bridge sites between two copper atoms. The copper atoms disorder intensely and some of them move onto the surface [Fig. 8(c)]. COMB3 thus induces an amorphous surface “oxide” layer as can be seen in Figs. 8(a)–8(c). A larger number of copper atoms in this amorphous surface oxide layer are no longer in an atomic environment close to their usual fcc environment. Oxygen molecules do not penetrate into the bulk during oxidation. Simulations at 300 and 600 K provide the same results with no molecular dissociation being observed. By contrast, with ReaxFF potential, oxygen molecules are mainly chemisorbed in hollow sites [Fig. 8(f)], and they do not disturb much of the positions of the surface copper atoms even at 600 K. Copper surface atoms are still in their fcc environment [Figs. 8(d)–8(f)].

We believe that the prediction of the ReaxFF potential is more realistic, because experimentally, some dissociation of the oxygen molecules is observed and no amorphous thin surface oxide layer is observed. However, we will also continue to use the COMB3 potential to see if an amorphous oxide layer affects surface slip irreversibility.

Neither CuO nor Cu₂O is nucleated in our simulations. The main reason is probably because only few oxygen molecules get dissociated during the currently available MD simulation times. CuO is monoclinic and Cu₂O is of the Strukturbericht C3 type (space group $Pn\bar{3}m$ ), with oxygen atoms occupying one tetrahedral site over four with respect to the copper fcc sublattice. The formation of CuO or Cu₂O out of fcc Cu and O thus seems more difficult than the formation of NiO out of Ni and O. The simulation times in this study are very short. Let us note that Cu₂O islands have been observed experimentally only at much larger Langmuir exposures and/or higher temperatures. ^{Reference Dubois57,Reference Zhou and Yang58}

2. Surface slip irreversibility of copper in oxygen environment

Copper surface slip irreversibility is first tested with COMB3 potential. As already mentioned, there is a great disordering of the copper atoms in the oxide layer with this potential. When a dislocation glides toward such an amorphous oxide surface, a surface step cannot be created in the same way as in the inert environment case or in the case of the nickel–oxygen system where nickel atoms in the oxide layer essentially maintain a crystalline fcc structure (see Sec. III and IV with a. 5, 6 and 7). What one now gets is shown in Figs. 9(a) and 9(b). We can propose the following interpretation. Although the concept of dislocation was developed by Vito Volterra in 1905 in the context of continuous elastic media such as gelatin or synthetic rubber, ^{Reference Volterra59} a dislocation in a crystal is a localized defect defined as a structural defect with respect to a perfect atomic crystalline order (regular ordering of atoms). Hence, in our case, when the edge 〈110〉/2 dislocation “reaches” the amorphous surface layer, the two extra atomic half planes which “correspond” to that dislocation cease to transfer their crystalline character and simply push outward the atoms of the amorphous layer. The extremity of the surface step which is correspondingly created simultaneously disorders and joins the amorphous layer, as shown in Fig. 9(b). As an obvious consequence, these dislocations cannot glide back under opposite loading and are fully irreversible.

FIG. 9.

Three snapshots taken shortly (about 50 ps) after: (a) the glide of two similar dislocations toward a disordered surface oxide layer at 300 K with COMB3, (b) a third dislocation has been added, and (c) three dislocations of opposite signs have been inserted in the fifth atomic plane above (as in Fig. 7). (c) The initiation of a nano-notch.

Even if we consider the oxidation behavior predicted by COMB3 potential as not too reliable with respect to real copper oxidation, it serves us to show that even in the special case of the creation of a soft amorphous layer which can be pushed outward by incoming dislocations, surface steps are initiated by the gliding of dislocations and are irreversible under opposite loading. A micro-notch can also be created under cyclic loading by introducing opposite sign dislocations, as can be seen in Fig. 9(c).

Surface step irreversibility and notch formation is also tested with ReaxFF potential. Inserted dislocations can get out of the surface oxide layer easily, because the surface copper atoms still retain more or less their fcc environment, which allows dislocations to transfer their character. This is similar with tests in inert environment and with tests with nickel and oxygen. The surface step will reconstruct and thus get reversible under opposite loading. If opposite sign dislocations are inserted continuously, a surface micro-notch forms.

C. Estimation of the total irreversibility

The atomistic simulations carried out in this work show a strong irreversibility of surface slips. Oxygen environment seems to have no major influence on the slip irreversibility and surface roughness compared with the case in inert environment. This is in agreement with the observations that the development of PSB surface relief is the same for fatigue in inert environment and in air. ^{Reference Thompson, Wadsworth and Louat7,Reference Greenfield20–Reference Witmer, Farrington and Laird24,Reference Basinski and Basinski60} The estimation about the surface irreversibility on the surface in inert environment that we have proposed in the previous work ^{Reference Fan, Hardouin Duparc and Sauzay5} should still be valid in oxygen environment. Here, we do a more complete estimation about the total irreversibility by coupling bulk and surface mechanisms. Both the Differt, Essmann, Gösele, and Mughrabi (DEGM) model ^{Reference Differt, Esmann and Mughrabi2} and the Lépinoux and Kubin (LK) model ^{Reference Lépinoux and Kubin3} are considered. The LK model is an improvement of the DEGM model.

As mentioned in Sec. I, the total irreversibility factor p is defined as the ratio between the irreversible slip and total plastic slip per cycle. It can be expressed using a weighted rule of mixture considering both edge and screw slip irreversibility factors:

$$p = {C_{\rm{e}}}{p_{{\rm{edge}}}} + {C_{\rm{s}}}{p_{{\rm{screw}}}}\quad ,$$

C _e and C _s correspond to the percentages of edge and screw slips contributing to the whole plastic slip, with C _e = 1 − C _s. The irreversibility factors p _edge and p _screw are defined as the total irreversibility factor, expect that those ratios are computed using only pure edge and screw dislocation slips, respectively. Taking into account both bulk and surface mechanisms, one has p _edge = p _edge (bulk mechanisms) + (1 − p _edge (bulk mechanisms)) × p _edge (surface mechanisms). And a similar equation holds for p _screw.

For edge dislocation slip, the irreversibility factor predicted by the DEGM and LK bulk models should be very small since the annihilation distance between edge dislocations are quite small. The estimation of “p _edge (surface mechanisms)” is based on our simulations. We observed that a surface step is reversible only when an opposite sign dislocation carrying out the same Burgers vector amplitude glides in a direct neighboring plane. A simple irreversibility calculation taking this neighborhood effect into account can be proposed: let us assume that a PSB lamella corresponds to a stack of N planes and that during one half-cycle n dislocations glide to that surface zone in a random way and create n steps. During a subsequent opposite half-cycle, one can assume that n opposite sign dislocations glide to the same surface zone (dislocation nucleation from the surface is not observed in our simulations). The numbers n and N are linked by the value of the plastic slip variation $\Delta {\gamma _{{\rm{pl}}}}{\rm{/}}2 = nb{\rm{/}}Nd$ with d being the interplanar (111) distance ( $d = a{\rm{/}}\sqrt 3$ , versus $b = a{\rm{/}}\sqrt {\rm{2}}$ ) and Nd varies from 0.4 to 1 µm in copper at room temperature. ^{Reference Basinski and Basinski60} The plastic slip amplitude $\Delta {\gamma _{{\rm{pl}}}}{\rm{/}}2$ is experimentally found to be of the order of 1%. ^{Reference Differt, Esmann and Mughrabi2,Reference Mughrabi61} If the material is “wavy,” with double cross-slips, the opposite sign dislocations during the opposite half-cycle will glide to the surface in a random way quite independent of the atomic glide planes of the first n dislocations. Using statistical combinatorics, a rough estimation of the surface slip irreversibility after one full cycle is thus p _edge (surface mechanisms) ≈ $C_{N - n}^n{\rm{/}}C_N^n$ , which gives a value between 0.5 and 0.75 for the copper measured data mentioned above. Knowing that the irreversibility factor estimated via bulk mechanisms considering edge dislocations alone is estimated to be very low, the value of p _edge considering both surface and bulk mechanisms can still be considered to be between 0.5 and 0.75.

Contrary to edge dislocations which disappear as such when they emerge at a free surface (one no longer has extra atomic planes in a bulk environment), screw dislocations always exist when they glide along surface where they create surface steps. A screw dislocation can thus easily glide back under reverse loading, even in case of reconstruction, because the surface step energy may be decreased by the gliding back of the screw dislocation. This means that surface steps created by screw dislocations are highly reversible [p _screw(surface mechanisms) ≈ 0] when there is no cross-slip. In wavy materials where cross-slip can occur easily, a screw dislocation can be easily annihilated by an opposite sign dislocation (by cross-slip). ^{Reference Smeenk, Tromp, Van Der Veen and Saris49} The irreversibility caused by this annihilation mechanism in the bulk is 0.3 in copper at room temperature following DEGM. ^{Reference Differt, Esmann and Mughrabi2} A larger value of 0.7 is found following the LK model, with similar considerations about the annihilation of screw dislocations by cross-slips and correlated movements of edge dislocations but by using a more refined model. As for screw dislocations, the complementary irreversibility induced by surface mechanisms is expected to be small, and the value of p _screw considering both surface and bulk mechanisms can still be considered to be about the bulk values, 0.3 (DEGM) or 0.7 (LK).

The total irreversibility factor amounts to p = C _e p _edge + C _s p _screw. For the sake of clarity, the different irreversibility factors are listed in Table II, depending on the edge/screw dislocation character and bulk/surface model.

TABLE II.

Synoptic table of the irreversibility factor calculations, from the DGEM and LK bulk models and from our MD computations for the special case of single crystal copper at room temperature. C _e and C _s are the fractions of slip carried out by the edge and screw dislocations, respectively. The total weighted irreversibility factor p is p = C _e p _edge + C _s p _screw.

It is found that our coupled total irreversibility factors considering surface mechanisms are two to three times higher than the values estimated by the DEGM bulk model or the related LK model. It should be noticed that the environment itself (oxygen molecules) is less influent than the pure free surface conditions investigated in detail. ^{Reference Fan, Hardouin Duparc and Sauzay5}

Similar calculations have been carried out for nickel using our irreversibility model considering both bulk and surface mechanisms, and oxygen environment. Similar values of p are found as those in copper. The estimated value obtained with the (LK + MD) model is not far from the value obtained from direct experimental measurement of surface steps by atomic-force microscopy on nickel cycle in air in the PSB regime, viz. 0.8. ^{Reference Weidner, Sauzay and Skrotzki4} The simulation and experiment conditions are consistent (material, temperature, plastic slip, oxygen environment), except that concerning the partial pressure values compensated by the difference in considered time scales.

Concerning planar materials where cross-slip is difficult, bulk annihilation of screw dislocations should be rare. This may decrease the irreversibility. Much lower irreversibility factors have been found in planar materials (see Table I in Basinski and Basinski ^{Reference Mughrabi61} and Mughrabi ^{Reference Mughrabi62} ). Further simulations need to be carried out and compared with experimental results obtained in such materials. In particular, additional MD computations should allow to check in more detail the behavior of surface steps induced by screw dislocation glide in both inert and oxygen environments.

Finally, our theoretical investigations of environmental effects on surface evolution, slip irreversibility, and nano-notch formation allow us to revisit the different models proposed in literature and described in the introduction section [mechanisms (i) to (v)]. Following our investigations on copper and nickel, oxygen adsorption and even dissociation do not prevent rewelding at all [mechanism (i)]. It has been proved that no emission of dislocation from free surfaces occurs on contrary to the proposal of Fine and co-workers. ^{Reference Sriram, Fine and Chung12} If a surface step is formed following the emergence of one-edge dislocation, then the step is erased during an opposite half-cycle if an opposite dislocation glides toward the free surface along a directly closed atomic plane. The surface oxygen atoms (adsorption + dissociation) or molecules (adsorption) do not hinder rewelding at all, as shown by the comparison of our new results in oxygen environment with the ones obtained in inert environment. ^{Reference Fan, Hardouin Duparc and Sauzay5} Therefore the hindering of rewelding largely discussed in literature ^{Reference Thompson, Wadsworth and Louat7,Reference Sriram, Fine and Chung12} [see mechanism (i) in the introduction section] does not seem to occur, at least for copper and nickel at room temperature. Interestingly, our simulations show that diffusion occurs quite quickly, and oxygen atoms penetrate until a few atom layers in the case of nickel at RT, as using ReaxFF potentials, without changing much the observed surface slip mechanism. We observe some effects of dislocation glide on oxygen penetration because cyclic glide seems to trigger oxygen penetration in agreement with the proposal of Fujita ^{Reference Fujita, Drucker and Gilman13} and the measurements of Fine and co-workers ^{Reference Sriram, Fine and Chung12} [model (ii)]. An increase of the friction stress is probably induced by oxygen penetration in the case of nickel, but once more that does not affect much surface relief evolution with respect to the evolution predicted in inert environment, ^{Reference Fan, Hardouin Duparc and Sauzay5} in contrary to the proposal of Fujita. ^{Reference Fujita, Drucker and Gilman13} An amorphous “oxide” subnano layer (copper, COMB3) does not affect the major trends concerning surface relief cyclic evolution. Further work is required to investigate the influence of oxygen adsorption/dissociation and eventually penetration on surface and fracture energies [mechanism (iii)] and nucleation of dislocation from notch tips [mechanism (iv)]. In agreement with models considering inert environment, ^{Reference Differt, Esmann and Mughrabi2,Reference Fan, Hardouin Duparc and Sauzay5} notches are produced due to slip irreversibility. Even if the notch formation is not much changed by oxygen environment, the transformation of notches into microcrack through cleavage [mechanism (iii)] or dislocation nucleation at notch tips may be environment dependent. Further MD simulations similar to the ones proposed in Fan et al. ^{Reference Fan, Hardouin Duparc and Sauzay5} and in this work may be useful to evaluate the effect of oxygen environment on mechanisms (iii) and (iv). The last model mentioned in the introduction, ^{Reference Shen, Podlaseck and Kramer17} corresponding to the interaction of pile-ups or slip bands with an environment-induced surface layer (v) would require the formation of a crystalline oxide layer to be efficient, which requires higher temperature/time than the ones investigated in our work.