Nanowire Stretching by Non-Equilibrium Molecular Dynamics

Non-equilibrium Molecular Dynamics (NEMD) simulations of a stretched Lennard-Jones (LJ) model single crystal nanowire with square cross-section are carried out. The microstructural and mechanical properties are examined as a function of strain and strain rate. The instantaneous Poisson ’ s ratio and Young ’ s modulus are shown to be strongly time (strain) dependent from the start of the pulling process. The structural transformation as a result of straining initially involves the (100) layers moving further apart and then slipping at ca. 45 o when the shear slip stress along that direction is about 1% of the shear modulus, which is typical of plastic deformation of noble gas solid crystals, and in accordance with Schmid ’ s law.


Introduction
The elastic properties of crystalline and polycrystalline materials have been the subject of practical interest over many centuries, and this field was placed on a firm theoretical footing by Navier, Cauchy, and Poisson. [1] The strength and failure mechanisms of solids under applied stress are important aspects of this topic. In recent years the remit of this field has broadened to include metallic nanowires which are used in miniaturized electronic and micro-machine (MEMS) devices [2] The mechanical and electrical properties of a given nanowire determine its suitability for particular applications. Large surface-to-volume ratios and other scale effects can cause the mechanical response to straining of materials of nano scale dimensions to be quite different to that exhibited by the same material with the same shape but of macroscopic size.
Under certain conditions, a wire a single atom thick can be produced at the final stage of the pulling process of a thin metallic wire. [3] The difference in behavior is, in part, because macroscopic wires are polycrystalline and deformation takes place predominantly by grain sliding at their boundaries. The inhomogeneity of the local structure leads to a heterogeneous stress and strain distribution in the wire. [4] In contrast, nanowires are, in the main, composed of single crystals, so that deformation mechanism under stress is not available to them. Another difference is that the strain rates can be many orders of magnitude larger in MEMS devices than those on the macroscopic scale, in part because the natural frequency scales as $ m À1=2 where m is the mass of the device element. The conditions experienced by the device therefore overlap to a much greater extent with atomic scale structural and dynamical characteristics.
The key mechanical properties of macroscopic wires under tension are Poisson's ratio, v, [5] and Young's modulus, E. Poisson's ratio is particularly informative as it is a direct and visible indicator of the nature of the deformation of the wire on stretching. Poisson's ratio has even been used to characterize the stretching deformation behavior of a single molecule, such as DNA. [6] It is usually derived from strain rates that are infinitesimal on a molecular scale. For a wire made of structurally isotropic material and of square cross-section, Poisson's ratio (PR) is the ratio of the transverse strain, ϵ xx or ϵ 33 in response to a strain, ϵ zz , along the length of the wire, where z is the pulling direction, and x and y are the transverse directions. Classical elasticity theory applied to infinitesimal linear strain gives, where G is the shear modulus and K is the bulk modulus, both measured in the limit of zero frequency. On a macroscopic timescale the elastic constants can be taken to be their isothermal equilibrium values. Young's modulus, E, is the constant in a generalization of Hooke's law for longitudinal deformation [7] E¼ σ zz =ϵ zz ¼ ÀP zz ϵ zz ; where σ zz is the wire axial or longitudinal stress and P zz is the corresponding pressure tensor element (note the sign difference). The classical theory of elasticity gives, [8] The classical linear strain theory gives the bounds, À1 < ν < 0:5, which goes from infinite shear modulus to rubbery solids. Although, for auxetic foams departures from Eq. (2) have been observed in experimental studies. [9,10] As a complement to experiment, nanowire mechanical properties can be conveniently explored using Nonequilibrium Molecular Dynamics (NEMD) simulations, which give a full atomistic level description of the system. This methodology can be employed to explore the dependence of the response as a function of external conditions such as temperature and strain rate. Here, classical dynamics and a pairwise additive pair potential function, ϕ r ð Þ, as a function of the separation, r, between the atoms, are used. A number of NEMD calculations of nanowires have been carried out using model metal force fields which have a many-body component, for example the embedded atom model (EAM). Wire stretching and nano-welding NEMD simulations have been carried out on a number of model metals, including copper [11][12][13] nickel [14][15][16] silver, [17] gold, [18,16] and Au-Co and Pt-Co alloy nanowires. [19,20] Also NEMD simulations have been carried out on βÀSiC. [21] These previous studies have focussed on the breaking mechanism of the wire on stretching through to the "necking" stage just before the wire breaks. In these calculations the model wire was attached to two cross pieces at the two ends that are pulled apart, to model a real apparatus. Substantial concerted restructuring of the wire through a sequence of potential energy minima was observed during the pulling process.
Other NEMD simulations of model solid wires under tension have been carried out using periodic boundary conditions in the axial direction, [22,23] and a comprehensive review of the nucleation of plastic deformation in nanowires is given in Ref. [24].
The stretching of liquids by NEMD has been carried out previously. [25][26][27] Techniques that are designed to exploit the flow behavior of the liquid state were employed in those cases, which would not be applicable for the present solid systems. The objective of the present study is to characterize the initial stages of the response of the solid nanowire to stretching, which could eventually lead to the improved design of mechanical elements in small devices. Inter alia the applicability and relevance of classical elasticity theory to nanowires is explored. The Lennard-Jones (LJ) potential, ϕ r ð Þ, was used to represent the interaction between the atoms, where ϵ and σ are the characteristic energy of interaction and diameter, respectively. The separation between the centers of two of the atoms is r. This is a simple generic atom or molecule pair function which includes short range repulsion and medium-tolong range van der Waals attraction. Its use helps to establish "universal" trends that are not dependent on a specific metal interaction law. Despite its simplicity the LJ potential has been used to model metals by molecular simulation in the past, and there are LJ ϵ and σ parameters for metals in the literature. [28] The LJ potential has been used to model metals in molecular simulation many times (see Refs. [29][30][31][32]). This work concentrates on details of the NEMD computer method developed for this project, and in revealing the microstructural processes responsible for the elastic-to-plastic deformation in terms of the classical parameters of material deformation under strain. The interest is in discovering the extent to which plastic (i.e., non-recoverable) deformation in a truly nanosized single crystal resembles that of more macroscopically sized single crystals for which the major source is motion of dislocations through the crystal lattice. This takes place parallel to the highly packed crystal planes of the original lattice, which are known as "slip planes". The dislocations transport the strain through the crystal by a concerted action of local atom exchanges with low activation energy. The initiating or yield shear stress, τ c , acting along the slip plane depends on the relative orientation of the crystal lattice relative to the pulling direction, and it is the pulling tensile stress, σ c resolved along the slip direction on the slip plane that must exceed a certain value to initiate plastic deformation, which is quantified by Schmid's law. This states that, where ϕ is the angle between the normal to the slip plane and the direction of initiating tensile stress (z here), and λ is the angle between the slip direction and z. [33][34][35] The largest value of cos ϕcosλ (the "Schmid factor") for a cubic crystal is 0.5, where both angles are 45 o . The extent of slip and increase in length of a single crystal depends on the magnitude of the loading stress, the crystallographic structure, and the number of active slip planes in the shear stress direction. When a single crystal is deformed freely in uniaxial tension without boundary constraint, the gliding planes retain their original orientation. Whereas, with boundary constraints (the case modeled in this study) the planes rotate toward the tensile axis to facilitate sliding and retain the connected integrity of the crystal. In addition, the lattice can rotate in the slip plane by a process known as dislocation "twinning", giving the extended crystal specimen a "twisted" appearance when viewed from the side or along the axial direction. In Section 2, aspects of elastic moduli of bulk crystals are covered. In Section 3 the NEMD simulation methodology for wire stretching is described, and in Section 4 the simulation results of model wire stretching are presented and discussed. Conclusions are presented in Section 5.

Elastic Constants Theory and Simulations for Bulk Systems
In this section a comparison is made between the elastic moduli and viscoelastic behavior of model liquids and solids. The stretching of liquid filaments is a widespread activity in industry, for example, in polymer melt processing, and it is of interest for future work to make a comparison between the two. Also it provides a context in which to introduce the formulas used to calculate the infinite frequency elastic moduli of the solid state, which are of immediate importance for the present work.
www.advancedsciencenews.com www.pss-b.com A number of mechanical quantities of use in characterising a nonequilibrium or finite strain (strain rate) process can be obtained directly from equilibrium state molecular simulations. The elastic properties of fluids [36] in the infinite frequency and low strain rate limit are so-called "static properties" as they can be calculated from equilibrium structural features at the level of the pair radial distribution function, g r ð Þ. For an equilibrium (i.e., unstrained) liquid with a spherically symmetric g r ð Þ, the infinite frequency shear rigidity modulus, G 1 is defined as [37,38] where k B is Boltzmann's constant, and ϕ 0 dϕ=dr and ϕ 00 d 2 ϕ=dr 2 . The first term on the right hand side of Eq. (7) is the kinetic part and the second term is the interaction or configurational part. Similarly for the infinite frequency compressional modulus, K 1 , is [38] These expressions can be reexpressed in a form that is suitable for implementation in simulation. Equations (7) and (8) are then, where V is the volume containing N molecules, r ij is the pair separation between molecules i and j, ϕ 0 ij dϕ r ð Þ=dr ð Þ r¼r ij and ϕ 00 ij d 2 ϕ r ð Þ=dr 2 À Á r¼r ij . Also, h. . .i denotes a simulation average, which is a time average for Molecular Dynamics (MD) simulation. Equation (9) is the Born term. [39] Similarly for K 1 , Alternatively, for the LJ potential, [40] these elastic moduli can be written in terms of the average repulsive and attractive components of the potential energy of the system expressed as an energy per molecule, and where u r ¼ U r =N and u a ¼ U a =N, and U r is the r À12 part of the potential energy of the N molecules, and U a is the same quantity for the attractive, r À6 part of the potential energy. [40] Let P xy be the shear (or an off-diagonal) component of the pressure tensor, which is calculated at each time step in the simulation using, where m i is the mass of molecule i, v xi and v yi are the x-and y-components of its velocity, respectively, and r xij is the x-component of the separation between molecules i and j.
The mechanical definition of the zero-frequency shear modulus, G 0 is, where P xy hP xy i is the time average (in MD) of P xy , in response to an applied shear strain, γ. Also from statistical mechanics [41,39] which is a combination of a nonfluctuating (Born) term, G 1 and a fluctuating term. For an equilibrium fluid, h P xy i ¼ 0, G 0 ¼ 0 which provides the definition of G 1 in terms of the shear stress fluctuations from Eq. (15), and for a solid G 1 > G 0 > 0. The elastic constants of the FCC Lennard-Jones solid have been computed numerically. [42][43][44] The differences between the viscoelastic features of a liquid and solid bulk system are compared in Figure 1, which shows the shear stress relaxation function, C s t ð Þ for a bulk equilibrium The state points are coexisting on the melting line at T ¼ 1:00, [40] and N ¼ 2048. The liquid density is 0:9157 and the solid density is 1:0030. Constant total energy or NVE dynamics were used.
www.advancedsciencenews.com www.pss-b.com system, in the linear response or infinitesimal shear strain limit. The time-dependent function, C s t ð Þ, was calculated (and is the same as) the shear stress time correlation function, which can be proved from Green-Kubo theory, [45,46] Figure 1 compares C s t ð Þ for liquid and solid coexisting states along the melting line, noting that C s 0 ð Þ G 1 for the liquid. [47] The melting line was chosen as the temperature and pressure are the same for the two co-existing phases, and only the densities are slightly different. This provides two closely related coexisiting thermodynamic states, but with different symmetry. The simulations were carried out using the leapfrog form of the particle update Verlet algorithm to generate the evolution of the system through time and space. [46] The reported results are in LJ reduced units (i.e., ϵ for energy, σ for distance and the mass of the molecule, m). The LJ interaction was truncated at r ¼ 2:5. [46] The fluctuations in the shear stress are larger in the liquid than the solid for the closely related state points. The figure also shows that the time scales for decay of the shear stress fluctuations to zero are quite similar.
The zero frequency isothermal bulk modulus, K 0 is a thermodynamic property defined as, where ρ is the mean particle number density. The values of K 0 from equilibrium simulations at the state point, T ¼ 1:0 and ρ ¼ 1:1 was 87:91 in LJ reduced units of ϵσ À3 . An NEMD simulation using Lees-Edwards periodic boundary conditions, [48] at a fixed shear strain of γ ¼ 0:02 was carried out to determine G 0 at this state point, giving a value of 53.72. Using G 0 for G and K 0 for K in Eq. (2) for v and Eq. (4) for E gives ν ¼ 0:246 and E ¼ 133:8. The value for ν ¼ 0:25 obtained agrees with that expected from materials whose atoms interact via central potentials which are isotropically distributed in space or obey the Cauchy relations. To an extent this aspect of the present work follows on from that of Ho et al. [49] who investigated negative Poisson's ratio in bulk cubic materials along its principal directions.
The elastic properties of bulk liquids and solids are to a large extent quite well understood. The situation is less clear for single crystal nanowires of nanometer dimensions, aspects of which are explored in the next two sections.

Nanowire Simulation Methodology
The NEMD simulations mimicked the pulling of a nanowire composed of a face-centered cubic (FCC) single crystal of square cross-section with the long axis along the [100] direction. The system was forced increasingly away from the thermodynamic and mechanical equilibrium state toward failure. The model nanowire had periodic boundary conditions applied in the axial direction, and was free to relax in the directions orthogonal to the pulling direction as there were no periodic boundaries in those directions.
The simulation cell had the dimensions, S Â S Â n z S, in the x, y, and z directions, respectively. For most of the simulations there were 1728 molecules in the simulation cell, which was made from two cubes (i.e., n z ¼ 2) of 864 molecules combined in the z direction. There were 24 xy layers of atoms in the zdirection. A simulation was also carried out with n z ¼ 3 and 2592 LJ atoms in the simulation cell. The simulation strategy leading up to nanowire failure was carried out in three stages.
The first stage (I) of the simulation was to allow the bulk system to come to equilibrium. The second stage (II) was to remove the periodic boundary conditions in the x and y directions. This was accompanied by a randomisation of the molecule velocities every 30 time steps, to help damp out the effects of the transition in physical state between the stage I to II conditions. This procedure created a nanowire of square crosssection and effectively of infinite length, as the system was still periodic in the z direction. The wire was allowed to come to equilibrium during this stage. In both stages I and II the model system was not stretched. In the third stage (III) the model wire was stretched by increasing the periodicity length in the z-direction. The stage III configuration is illustrated schematically in Figure 2.
The initial periodicity length, S z 0 ð Þ ¼ n z S, was increased linearly with time, t, through intermediate values, The typical strain rate was 6:3 Â 10 À4 reduced units which corresponds to $ 10 8 s À1 for solid argon, which is many orders of magnitude greater than is typical of macroscopic fiber stretching experiments. The bottom boundary of the periodic cell was fixed at, z ¼ 0, while the upper boundary was gradually increased by a small incre ment each time step until a maximum extension, ϵ zz;m ¼ S z t m ð Þ À S z 0 ð Þ ð Þ =S z 0 ð Þ was achieved, and where t m is the maximum time of the stage III segment. This procedure was implemented to mimic as closely as possible the stretching of a real wire. The atom coordinates in the z-direction were not uniformly scaled, which would have been unrealistic. Density and other system property inhomogeneities should be allowed to evolve naturally without interference from the stretching mechanism. The implemented mechanism of stretching followed closely the experimental procedure, in which the wire is clamped at two points and pulled apart from there.
The system's temperature in stages I and II was controlled with the Nos e-Hoover thermostat, [50][51][52] with a time constant of three reduced units. During the stage III stretch of typically 75 reduced time units the system was not thermostatted but allowed to reach a time-dependent temperature in response to the stretching process, which would be a natural feature of the pulling process in the real world. Stages I-III were repeated N s times to improve the statistics of the calculated properties. After each stage III, the molecule positions and velocities saved at the end of the previous stage I phase were recovered as the starting point of the next stage I time period (of ca. 25 reduced time units), which was used to create a new starting state independent of the previous one for the stretching event. Therefore successive wire equilibration and stretching events commenced from statistically independent starting states. The number of such independent stretching trajectories, N s , was typically % 200. The density was 1.0 and the temperature was equal to 0.1 to ensure www.advancedsciencenews.com www.pss-b.com the starting system was in the crystalline part of the phase diagram. [53] The time step at the target temperature of 0.1 was 0.016 LJ reduced units. A typical simulation consisted of 5000 time steps in stage I, 25 000 time steps in stage II and 10 000 time steps in stage III. The velocities were randomized by reassigning them from uniformly distributed random numbers and then scaling them to the target temperature. A Maxwell-Boltzmann distribution of velocities was achieved after about 100 time steps.
Also in stage III the angular momentum of the wire around the axis was set to zero to prevent the wire from spinning around the z-axis, which it will because of the residual velocities and that bulk MD does not conserve angular momentum. This is unwanted for statistical mechanical reasons and also on a more practical level would have made the calculation of the lateral strains and Poisson's Ratio problematic. First the angular velocity ω, around the z-direction and through the center of mass was calculated from, where r i is the xy plane projection of the coordinates of atom i relative to the center of mass, and v i is the component of the velocity in the xy-plane of particle i. The correction to the individual x and y coordinate velocities was,

Nanowire Modeling Results
The atomic-scale structure of the model nanowire at various stages in the pulling process implemented by NEMD is presented in this section. At the temperature used, T ¼ 0:1, the wire is stable, and under a slight tension (at this temperature the original bulk system is under a negative pressure of À6:9 LJ reduced units). The radial distribution function of the bulk LJ solid with FCC crystalline form for the equilibrium state point T ¼ 0:1 and ρ ¼ 1:0 is presented in Figure 3, which shows the thermally broadened peaks characteristic of an FCC crystal. The reason why having a stage II in the simulation procedure is necessary is illustrated in Figure 4, which shows the pressure component, P zz as a function of time in stage III obtained from the Method of Planes (MOP) method for four xy-planes equally spaced in the z-direction. The simulation was carried out with a stage II of zero time duration, which means that the stage III data were accumulated immediately after the periodic boundary conditions in the x-and y-directions were removed. To obviate any masking by the stretching process itself, the strain rate and maximum strain over the 80 reduced times step stage III segment were given insignificantly small values of _ ϵ zz ¼ 1:3 Â 10 À10 and ϵ zz;m ¼ 10 À8 . Oscillations in the axial component of the pressure tensor across various xy planes along the wire are seen to be in phase and long-lived (rather like a  www.advancedsciencenews.com www.pss-b.com struck bell vibrating). This may be attributed to a damped oscillatory relaxation of the model nanowire cross-sectional area as a result of the "deletion" of the "boundary" pair interactions when the periodic boundary conditions were taken away. This phenomenon may share similar physics to that responsible for quench echoes in glasses and proteins, which describes the recovery of a system to an abrupt temperature change. [54][55][56] Figure 5 shows xz-plane projections of the wire structure at z-strain values, ϵ zz of 0.0, 0.025, 0.050, 0.075, and 0.1 (the maximum value). The strain rate used was _ ϵ zz ¼ 0:00063. The ϵ zz;m ¼ 0:1 or final time step yz-projection is also shown (in red on-line) at the far right hand side of the figure. The x and y lengths are compressed to enable six frames to be shown on a single figure. The layers initially stay more or less parallel and just move uniformly further apart during the pulling process. Then close to ϵ zz ¼ 0:1 the layers can be seen to no longer lie in the xy-plane but are reorientated to a certain extent toward the z-direction. There is also evidence of a "twist" in the wire in the two rightmost snapshots, which exhibit ordered and apparently disordered regions that "spiral" around the column. Note the absence of particles in the bottom of the frame at ϵ zz ¼ 0:075, which are visible in the top layer because of the periodic boundary conditions in the z-direction; they could be drawn at the bottom of the frame as well as the system is periodic in that direction.
One might expect the uncontrolled failure of the wire to appear when the nearest neighbor distance exceeds the point of inflection distance of the LJ potential. This critical strain might be viewed as being the maximum stress and recoverable strain that can be sustained before permanent deformation and yielding takes place. The nearest neighbor distance, r nn , at the ρ ¼ 1:0 density is 2 1=6 , and the strain ϵ zz required for r nn to equal the point of inflection of the potential at r nn ¼ 26=7 ð Þ 1=6 is 0.23 (from r nn ¼ 1:122 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ ϵ zz p ) which is much larger than is found necessary here for plastic deformation. This is because slip plane gliding is in fact the plastic deformation mechanism that occurs in single crystals, as discussed in the Introduction. Also as revealed in Figure 6 the slipping planes of atoms reorientate toward the tensile axis to facilitate sliding, and the lattice rotates giving it a "twisted" appearance, as predicted in the literature. [57,58] Figure 6 gives an example of the sideways view of the wire looking along the y and x directions, from the same system as in Figure 5. The scale is the same in all cartesian directions. The structural distortion takes the form of a pattern which shows the importance of the 45 o planes to the xy plane. Figure 7 shows the xy-projections of the atoms in four of the layers within the 36 (increasing numerically from bottom to top in the periodic cell) at the end of the stretching process, when γ ¼ ϵ zz;m ¼ 0:1 and for _ ϵ zz ¼ 0:00032. This simulation used n z ¼ 3, so the simulation cell dimensions were 1 Â 1 Â 3. The layers are quite close to the boundary near the bottom of the repeated simulation cell. The figure shows that the layers retain much of their original initial crystalline structure, but they show evidence of dislocations consisting of three closely spaced layers out of the plane of the image, which runs vertically through the xy layers of the wire near the repeat unit boundary.
The instantaneous strain was calculated by taking the average distance between the x-coordinates of the left and right boundary atoms for each layer, or, hs x i and then using, where s x 0 ð Þ is the corresponding quantity at the beginning of stage III. The overall strains orthogonal to the pulling direction, ϵ xx and ϵ yy as a function of ϵ zz were taken as the average of those of each of the layers. This procedure removes any systematic errors which may be introduced by instantaneous variations in the wire thickness along its length at t ¼ 0. The layer-averaged lateral strains, ϵ xx and ϵ yy and their average are shown as a function of the longitudinal strain, ϵ zz in Figure 8 for _ ϵ zz ¼ 0:00063 and _ ϵ zz ¼ 0:00032, using the average response of 200 statistically independent stage III segment trajectories. The figure shows  that the wire becomes narrower as it is pulled, as the lateral strains are negative. The ϵ xx and ϵ yy are also statistically indistinguishable, at least up to a strain of ca. 0.08. Then there is a sharp increase in the magnitude of the strain which corresponds to the point where the layers start to reorient noticeably toward the pulling direction. The wire becomes thinner for that reason. There is some statistical variation between the behavior of the x and y direction in this strain region, which reflects a breakdown of the original structure. Data for the average lateral strain from a simulation carried out with half the strain rate is also shown on the figure, which indicates that there is no statistical difference between the two responses. This suggest that the response is quasistatic for the strain rate regime considered here. The sharp increase in the lateral strain after about ϵ zz ' 0:08 is because the wire becomes thinner as a result of the reorientation to an extent in the z direction. The reason the layer reorientation takes place is presumably because this is a means by which it can retain its connectedness at high levels of extension.
As the x and y strains are equivalent by symmetry, the Poisson's ratio is conveniently defined here as, The strain dependent Poisson's ratio derived from the data in Figure 8 is shown in Figure 9. It starts off with a value of ca. 0.20 and increases to ν ' 0:29 at ϵ zz ' 0:02 before decreasing to a value of ca. 0.23 at ϵ zz ' 0:07. From ϵ zz ' 0:08 to 0.10 the PR increases sharply which corresponds to a more rapid rate of increase of the axial strain as already pointed out. The PR may be a useful order parameter to describe the nonequilibrium potential energy landscape of the wire, as may be seen in this figure. Figure 8 also shows that Poisson's ratio is statistically independent of strain rate in the range considered, strongly suggesting, again, that the system is in the quasistatic limit, at least in the early stages of the stretching process.
The instantaneous diagonal elements of the pressure tensor were computed by the Method of Planes route. [59,60] The P zz component was obtained at the boundary (i.e., z ¼ 0 or z ¼ S z   www.advancedsciencenews.com www.pss-b.com by virtue of the periodic boundary conditions), and for three other planes inside the cell at z ¼ S z =4, z ¼ S z =2, and z ¼ 3S z =5. The P xx component of the pressure tensor was computed from interactions crossing the yz plane passing through the middle of the wire. These pressure values are shown as a function of ϵ xx t ð Þ for _ ϵ xx ¼ 0:00063 in Figure 10. The P xx e xx ð Þ is seen to be statistically zero as the model wire is surrounded by a vacuum. The z-component of the pressure tensor starts off negative (i.e., it is slightly under tension to begin with). The lateral strain is initially linear with axial strain up to ca. ϵ zz ' 0:02. Softening is evident for larger axial strains. Up to ϵ xx ' 0:04 all planes through the wire give statistically the same value of the pressure, but beyond that point P zz measured at the periodic boundary goes more negative more rapidly with axial strain. This could be viewed as boundary region hardening, presumably because the straining is applied at the periodic boundary, and it is in this region where the wire will at least for a period develop a structure which resists further stretching. The boundary region undergoes work hardening while the central region of the wire is still undergoing plastic deformation softening. It can be seen that for ϵ zz > 0:08 the system as a whole softens, evident in a decrease in the magnitude of P zz . The top and bottom layers undergo different localized structural changes from the rest of the system, resulting in an inhomogeneous distribution of stress. At this time (or equivalently global z-strain) the wire stretches to different extents along its length. The largest negative stresses are at the periodic boundaries (open square symbols on the figure). This is evident in Figure 5, and this state is a precursor to much larger structural rearrangements which take place afterwards for global strains corresponding to 0.08 and higher. This is an intermediate metastable stage which provides the free volume for larger restructuring to take place. The nz ¼ 2 and 3 cases had 24 and 36 layers of atoms, respectively. They both showed the "twist" structure at a global strain of 0.1, but it took longer to develop fully in the 36 layer case. At a global z-strain of 0.1 the twist region was more localized near the periodic boundaries in the nz ¼ 3 case.
Strain hardening prior to material failure is often observed in metallurgical experiments, and one supposition is that regions of the crystal are rotated relative to the rest which act as a break on the slipping of the surrounding crystal, enabling the whole system to withstand a greater load than the perfect crystal. [34] The transition from linear to nonlinear tensile pressure which occurs at ϵ zz ' 0:02 corresponds to a decrease in pressure from the starting values of about 1.5 reduced units (see Figure 10). If we assume that ϕ and λ are both 45 o , then Schmid's formula gives a yield slip shear stress of ' 1:5=2 ¼ 0:75 reduced units, and as the shear modulus of the crystal is ca. 54 (see above), then the ratio of the yield stress to shear modulus is 1.4% which is typical of the value found experimentally for solid argon. [61,62] This is an upper limit as the Schmid factor can be smaller than 0.5 (e.g., 0.3) but the conclusion is not materially altered by changing the value of the Schmid factor by a physically reasonable amount (the angles ϕ; λ have not been determined in this study). Figure 11 shows the corresponding Young's modulus, E, as a function of strain computed using the formula in Eq. (3). The figure shows that E e zz ð Þdecreases initially with increasing strain from a very large initial value, which levels out at values ' 80 in Figure 9. Poisson's ratio v for the strain averaged over all layers for ϵ xx ("x" on the figure), ϵ yy ("y"), and the average of ϵ xx and ϵ yy values ("(x þ y)/2") (see Eq. (1)) for the applied conditions, _ γ ¼ 0:00063 and γ m ¼ 0:1. The data denoted by "[x þ y]/2" on the figure is the average of the ϵ xx and ϵ yy for _ ϵ zz ¼ 0:00032 and ϵ zz;m ¼ 0:1. Figure 10. Pressure along the axis in the z-direction and across the yz plane calculated by the Method of Planes (MOP) method across selected planes along the wire. [59,60] The P zz component was computed at the boundary (i.e., z ¼ 0 or z ¼ S z ) by virtue of the periodic boundary conditions, and for three other planes, z ¼ S z =4, z ¼ S z =2, and z ¼ 3S z =5, where _ ϵ zz ¼ 0:00063 and ϵ zz;m ¼ 0:1. The P xx component of the pressure tensor was computed from interactions crossing the yz plane passing through the middle of the wire. The data labeled ½P zz is for the lower strain rate of, _ ϵ zz ¼ 0:00032 and ϵ zz;m ¼ 0:1.
www.advancedsciencenews.com www.pss-b.com the strain range 0:02 À 0:04. This part of the profile does not depend on where the plane for the MOP pressure is chosen. Then above ϵ zz ' 0:04 the value of E at the boundary increases sharply, while those computed from planes within the periodic cell start to decrease. It is at the boundary between the periodic cells where the stretching process is imposed on the system, so one might expect these region to behave differently when the system has large strains. The trends are in close agreement for the two strain rates considered. This sequence of variable strain behavior is broadly consistent with the known behavior of strain hardening of FCC crystals, [34] in which initially there is facile gliding with low hardening, then increased hardening (here near the periodic cell's boundary), which is nearly independent strain rate. This is followed by significant material weakening, further hardening, and then ultimate failure over a narrow strain range above ϵ zz ' 0:1. Figure 12 shows a map of the dislocation or "defect" locations represented by large colored circles superimposed on a three dimensional view of the lattice of atoms which are shown as small black dots. The dislocation site is identified using a measure of the spatial anisotropy of the positions of the surrounding neighboring using Eq. (10) of Ref. [63], which is also given in the figure caption. It measures the extent of breaking from the inversion symmetry of nearby atoms around a particular atom found for a perfect lattice. Snapshots of the system atoms at zz strains 0.07 (left image) and 0.10 (right image) are presented. The defects first form near the top and bottom of the repeated cell, which is what might be expected because the cell pulling apart is implemented at the boundary. The local structural anisotropy there reflects the anisotropic repositioning in the z-direction of the atoms on the two sides of the border between the cell and its periodic image. At a larger strain of ϵ zz ¼ 0:10 the dislocation defects have spread, starting at the lower boundary, in two lines which move through the solid and eventually populate most of the cell.
During the stretching process the thermodynamic properties change as a function of strain. Figure 13 shows the strain Figure 11. Time or strain-dependent Young's modulus, E e zz ð Þ. Most of the data plotted is for a strain rate of _ ϵ zz ¼ 0:00063 and ϵ zz;m ¼ 0:1. The data labeled z ¼ SZ=2 ½ is for _ ϵ zz ¼ 0:00032 and ϵ zz;m ¼ 0:1.   The temperature is given, as is the kinetic energy per particle, the potential energy per particle and the total energy per particle. The temperature profiles for _ ϵ zz ¼ 0:00063 and 0:00032 are given. The temperature data for _ ϵ zz ¼ 0:00063 is multiplied by 100 on the figure, and referred to as, "Temperature Â 100". The corresponding plot for _ ϵ zz ¼ 0:00032 is indicated by "[Temperature Â 100]". www.advancedsciencenews.com www.pss-b.com dependence of the nanowire temperature, T, obtained from the kinetic energy per particle, the potential energy per particle, and the total energy per particle. The figure shows that the temperature steadily decreases with increasing strain until the region of incipient failure (i.e., ϵ zz > 0:08), when it increases sharply. The particles start to get further apart up to that region, so that for ϵ zz < 0:09 a steady decrease in the magnitude of the potential energy and total energy is observed. There is hardly any statistical difference between the temperature profiles for _ γ ¼0:00063 and _ γ ¼ 0:00032.

Conclusions
Non-equilibrium Molecular Dynamics (NEMD) was used to model the pulling apart of a model nanowire. As with some previous simulations designed to model specific metals, [22,23] the system was periodic in the axial direction. Some properties of infinite and finite periodic systems in the wire geometry can differ if the period in axial direction is not large enough. [64] Nevertheless, the advantage of this procedure is that finite size effects are minimized, and the infinite wire limit could be estimated by extrapolation of the results for different periodicity lengths to infinite periodicity length.
In the procedure adopted here, two sample preparation stages are undergone before inducing tensile stress. The first involves equilibration of a bulk MD crystal system, where the unit cell is a right parallelepiped, taking the long axis to be along the z-direction. Periodic boundary conditions (PBC) are applied in all three cartesian directions. Then the PBC are removed in the x-and y-directions, and the newly formed LJ wire shape is allowed to settle to a new equilibrium state. This is, for example, to allow the vibrations that are created in the model wire by removing the PBC to die away. At the end of the nanowire equilibration phase the angular velocity about the wire axis is set to zero. Stretching of the equilibrated wire along the z-direction is subsequently carried out by moving the top boundary of the (repeated) simulation cell at a constant velocity along the z-direction. The MD cell remained periodic in the z-direction, and the system was not thermostatted in this last stage. Distortion of the model wire took the form of the (100) layers orthogonal to the wire axis moving further apart initially. Plastic deformation took place along glide planes at 45 o to the 100 planes, and then in the final stages these planes partly reorientated toward the pulling direction to form a stack of displaced sliding layers. This allowed the wire to retain its connected integrity. The model wire also appears to have "twisted" about the wire axis (see the two most rightmost frames in Figures 5 and 6). When the applied strain approached a value of 0.1 the xy-plane layers had reorientated noticeably. This caused a decrease in the thickness of the wire, an increase in layer averaged Poisson's ratio, and an increase in Young's modulus. There is evidence of strain hardening where the deformation of the wire is introduced at the periodic boundaries. In general it might be expected that the plastic deformation and failure mechanism could be sensitive to strain rate and crystallographic orientation.
Poisson's ratio and Young's modulus continually varied during the pulling of the model wire. This indicates that the potential energy landscape of the wire continuously changes during stretching. Also it might be concluded from this study that a time dependent (e.g., oscillatory) strain rate could be used to control the mechanical response and effective elastic constants as a function of strain.
Although more extensive investigations using different input parameters are possible, the structural transformations observed appear to be broadly consistent with those discussed in the single crystal plastic deformation literature. A more detailed study would be required to discern any differences, which one might expect, as a consequence of the large surface-to-volume ratio of the model nanowires, and perhaps the larger strain rates that are achievable on the nanoscale.
Future studies could explore the effects of the width and length of the wire on the plastic deformation mechanism, as it tends to the bulk system limit. Also, by removing the periodic boundary condition in the z-direction, and replacing it by a movable solid constraint, it would be possible to explore further the origins of the wire "twisting" about the wire axis seen in Figure 5.