Car–Parrinello Molecular Dynamics Simulations of Tensile Tests on Si

ARTICLE pubs.acs.org/JPCC

Car
Parrinello Molecular Dynamics Simulations of Tensile Tests on SiÆ001æ Nanowires Al'ona Furmanchuk,† Olexandr Isayev,‡ Tandabany C. Dinadayalane,† and Jerzy Leszczynski*,† †

Interdisciplinary Center for Nanotoxicity, Department of Chemistry and Biochemistry, Jackson State University, Jackson, Mississippi 39217, United States ‡ Department of Chemistry, Case Western Reserve University, Cleveland, Ohio 44106, United States

bS Supporting Information ABSTRACT: Theoretical simulations of tensile tests on SiÆ001æ nanowires have been carried out using Car
Parrinello molecular dynamics. H-passivation was used to model experimentally occurring passivation in Si nanowires. First-principle molecular dynamics simulations at ambient temperature reveal the governing role of size, overall shape, and composition of the surface layer for the mechanical properties. Our results indicate that SiH2 groups in the outer layer and the octahedral shape of the wire soften Young’s modulus and allow wire to handle larger transverse strains than SiH groups in wires with the tetrahedral shape. The importance of the overall shape of the wire has been discussed by comparing the behavior of surface layers of {100} and {110} facets. The presence of the {100} facets helps to relax the transverse strain during tension. On the basis of changes in structural parameters, we have presented the schematic motion of Si atoms in core and surface layers before the fracture appeared.

1. INTRODUCTION A silicon nanowire (SiNW) is considered as a promising alternative for a variety of nanodevice applications1
3 due to their unique mechanical and other properties.4
11 One of the reasons for a good performance of nanowires is the extremely high value of their surface area to volume ratios. Since the mechanical properties of surface layers are significantly different from those in the bulk, they will dramatically change the overall performance of semiconducting devices. Numerous experimental and theoretical studies showed the size-dependent character of the mechanical properties that are evaluated by Young’s modulus. It has to be noted that the presence of surface stress,12 contaminations,13 oxidation,14 the way of applying load (resonance versus bending or versus stretching/compression),15 and fabrication-induced16 or equipment-induced defects17,18 may affect the experimental measurement of Young’s modulus. Even the minor appearance of any kind of listed effects may result in different values of Young’s modulus.19
23 An increase of the surface-to-volume ratio of the SiNW decreases the accuracy since handling of the sample becomes tough. Precise control over the strength of SiNWs may be gained by performing theoretical simulations of mechanical tests. The right choice of a theoretical approach should be made to achieve an acceptable compromise between the simulation time and correct description of the quantum-mechanical effects caused by the surface. Numerical approaches allow us to model systems of almost experimental length, but there is no quantum-chemical insight into atomistic interactions. Empirical Si potentials24,25 for r 2011 American Chemical Society

the simulations of large-size systems give only an implicit description of the electronic effects and are incapable of describing the surface effects properly. On the other hand, a semiempirical tightbinding approximation26 lacks in multicentered integral calculations, although it is much better in performance. A good alternative to the empirical and semiempirical approaches would be the utilization of first-principle density functional theory based quantummechanical Kohn
Sham equations.27
30 However, the models used in tight-binding- and first-principle-based simulations are slightly different from the procedures used in experiments. For example, in the simulation of a free-standing sample, which has been uniaxially stretched during the tensile test experiments, the force is applied at both ends of the finite size sample (or the size of the super cell is changed in periodic systems). However, subsequent geometry optimizations of the system do not account for temperature effects explicitly. Therefore, the resulting process of stretching deviates from the one performed at ambient conditions in the laboratory. To explicitly include temperature effects, the density functional based, first-principle Car
Parrinello molecular dynamics (CPMD) simulations may be employed. These firstprinciple-based simulations properly treat the quantum chemical effects caused by atomic interactions at the surface layer. The Car
Parrinello scheme has been successfully used in studies of molecular rupture.31
36 Received: February 28, 2011 Revised: May 6, 2011 Published: June 07, 2011 12283

dx.doi.org/10.1021/jp201948g | J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C Similar to the theoretical challenges, there are experimental difficulties in growing SiNWs as described in the literature.19,37 To our knowledge, only two experiments, which were performed by Kizuka et al.19 and by Collins et al.,37 are available for a comparison with the theoretical data of the Æ001æ SiNWs. However, none of them has been performed strictly following the procedure used in the tensile test. In the first experiment, hardly defined shape and partial passivation of the sample do not allow for a direct comparison with the theoretical results. The large size of a ball-shaped specimen, which was used in the second experiment, is different from the model sizes studied theoretically. The mechanical properties of such specimen are rather comparable with bulk Si. Importantly, the values of Young’s modulus from the firstprinciple studies, which were reported earlier, on the Æ001æ SiNWs are not consistent. The only common feature in all those results is a general trend. The reason for the discrepancies is the utilization of different assumptions for the calculations of geometrical parameters. For example, taking a mean diameter in place of the wire width or a different definition of the cross sectional surface area28,30 may result in a 0.1
0.2 nm difference in the wire width and up to 23
24% reduction of the value of Young’s modulus. Even taking into account the presence of some discrepancies within the theoretical studies, we notice significant deviation between the experimental and any ab initio theoretical data. Such a difference may be attributed to the difficulty in modeling the effect of local surface strains due to uneven contours of the nanowires. Partially, it may be due to the H-passivation used for dangling bonds in simulations. However, the H-passivation is generally assumed as the simplest approach to the oxide layer passivation of the experimentally grown SiNWs. It seems to be present in the wet and chemically etched wires that were not oxidized.5,16,38
41 The effect of H layer on the values of Young’s modulus of the Æ001æ SiNWs under tension should be investigated to gain knowledge of the mechanical properties of oxide-free Æ001æ SiNWs. The present work is aimed at examining how shape and size of the Æ001æ SiNW affect its mechanical properties. The Car
Parrinello molecular dynamics simulations of the tensile tests were performed for the SiNWs of different widths and shapes. Our results are compared with the earlier first-principle study.28 New features of the size-dependent trend of the mechanical characteristics in small SiNWs have been revealed. Differences in the mechanisms of breaking, which are specific to the shape and size of SiNWs, are also discussed. The present simulations are important because of (i) the lack of experimental data for the tensile tests involving Æ001æ SiNWs and (ii) the poor agreement between the available experimental and theoretical data.

2. METHODS In accordance with experimental observations8,38,42 and Wulff construction,43 we model six nanowires by cutting them from bulk silicon using the low free-energy {100} and/or {110} facets. The longitudinal axis was oriented in [001] direction in all cases. The dangling bonds at the surface were passivated by H atoms in the following way: monohydride type of bonding formed over the {110} facets, while dihydride configurations formed over the {100} facets.44 No complex of SiH3 is observed over any facet of the wires. All simulations were performed using the standard Car
Parrinello molecular dynamics scheme.45 The Becke, Lee, Yang, and Parr (BLYP) exchange correlation functional with

ARTICLE

Figure 1. Cross sections of equilibrated Si nanowires. Typical direction of {110} and {100} facets labeled on sample B2 as an example. Color of Si atoms is yellow, and H atoms are given in gray. For the sake of clearness, the number of atoms within the cross section is as follows: A1, Si9H12; A2, Si25H20; A3, Si49H28; B1, Si21H20; B2, Si57H36; B3, Si109H52.

the generalized gradient approximation was used.46,47 For all atoms, we used the nonlocal norm-conserving Troullier
Martins pseudopotentials to describe the valence
core interactions.48 All pseudopotentials have been transformed to a fully nonlocal form using the scheme of Kleinman and Bylander.49 The valence electronic wave functions were expanded in a plane wave basis set up to an energy cutoff of 50 Ry. Hydrogen nuclei were treated as classical particles with the mass of the deuterium isotope.50 Simulations were performed at constant volume using a fictitious electronic mass of 400 au.45 The equations of motion were integrated with a time step of 5 au (0.121 fs). Each passivated wire was placed in an orthogonal cell with the cell walls separated from the wire by 10 and ∼5 Å correspondingly in longitudinal and transverse directions. All atomic positions were relaxed for 1 ps to reach the equilibrium structure. During this run, the temperature was controlled by uniformly scaling the atomic velocities to the target simulation temperature of 300 K whenever the tolerance interval of 50 K was exceeded. The system was then equilibrated at 300 K for about 1 ps using a Nos
e-Hoover chain thermostat51 with a chain length of four and frequency of 600 cm
1. In this work, six different widths of SiNWs are considered. The corresponding cross sections of the equilibrated wires are depicted in Figure 1, and the structural parameters are listed in Table 1. The rupture simulations were carried out in the NVT ensemble using the same Nos
e-Hoover chain thermostat condition. The “Blue Moon” ensemble method52
59 was adopted to stretch the Si nanowires. Thus, constrained dynamics was performed by adding a constraint to the Car
Parrinello Lagrangian LCP LCP ¼ LCP þ λξ ½ξðfRi
Rj gÞ
ξ0 

ð1Þ

where λξ is the Lagrange multiplier relevant to the constraint ξ({Ri
Rj}) that fixes the reaction path or, in our case, the distances between the ith and jth particles. Several pairs of the 12284

dx.doi.org/10.1021/jp201948g |J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C

ARTICLE

Table 1. Type and Number of Atoms Comprising SiNWs, Wire Width (nm), Initial Averaged Length (L0, nm), and Surface to Volume Ratio (nm
1) Calculated from Equilibrated Samplesa wire code

# of atoms (Si:H) within L0

wire width after equilibration, nm

L0, nm

surface/ volume, nm
1

A1

40:56

0.662

2.258

6.23

A2

109:92

1.041

2.172

3.82

A3

212:128

1.435

2.265

2.79

B1

90:90

1.068

2.172

3.75

B2

244:162

1.612

2.172

2.38

B3

466:234

2.164

2.257

1.77

a

See Figure 1 for structuresb and wire codes. b A1, A2, and A3 have a tetrahedral shape of the cross-sectional area; B1, B2, and B3 have an octahedral shape of the cross-sectional area.

second row silicon atoms from both ends of the wire were chosen as ith and jth particles. The ξ0 is
5.82  10
15 m fixed during the CPMD run. Such a stretching simulates how a tensile stress is applied in the tensile test experiments.19,60,61 The pulling velocity in our simulations (∼ 48 m/s) is faster than that typically used in experiments. However, it is slow compared to the vibrational movements of the nuclei. The total force for all velocitydependent constraints was collected during each step of the rupture run and was used for the estimation of the stiffness of SiNW under external tension. The tensile stress
strain curves are plotted for all Si nanowires that we simulated. The stress, σ, is calculated using the applied force divided by the original cross sectional area, which can be calculated by averaging the area of the outermost hydrogen atoms of the Si nanowire. The reason for using hydrogen atoms was well explained by Lee and Rudd.28 The strain, ε, is the elongation of SiNW divided by its initial length. Following Hooke’s law, Young’s modulus, E, is calculated from the slope of the initial linear region of each stress
strain curve by applying a linear fitting to the original data. Owing to high fluctuations on stress
strain curves, we applied a Savitzky
Golay filter62 to smooth the data. Thus, we are able to visualize the ultimate tensile and breaking strains and the ultimate tensile strength. The best improvement in data variance was reached when the frame size was set to 8000 points and the polynomial degree was equal to 1. To show how interatomic bonds realign during the stretching, we also calculate Poisson’s ratio, υ, for the SiNWs considered in this study. We define Poisson’s ratio in the elastic linear region as given below υ¼


ðd
d0 Þ εd0

ð2Þ

where d0 and (d
d0) are the original width of SiNW and its change during stretching, respectively. All geometrical parameters presented in this article are measured within the L0 length of the SiNWs and averaged for each step of the simulation. We would also like to mention the limitations that one may face during similar simulations. Our approach would not be quite reliable for the wires whose widths are smaller than 0.5 nm. Even for the 0.66 nm wire, we can still experience the impact of twisting. Due to the way the constraints are applied, one has to work with the specimens that undergo less twisting. However,

mechanical properties can be reliably determined in the region of lower strains for the Si nanowires whose widths are equal to at least 1 nm.

3. RESULTS AND DISCUSSION 3.1. Zero-Strain Structures. We would like to start with the comparison of key structural characteristics of our results of the equilibrated systems with the previously reported ones.27,28,39,63
65 For example, the correct description of dihydride configurations of the H-passivation would be a good test for the methodology justification. From Figure 1, one may see that Si-dihydrides have a rather canted character, which means a large enough separation between the H atoms of the two nearest SiH2 groups. Such a separation reduces the repulsion between the H atoms of the two nearest SiH2 groups and thus stabilizes the structure. The preferential stability of canted SiH2 over symmetrical ones was well established by Northrup for the bulk Si surfaces.63 It was also supported by thorough theoretical27,28,64,65 and experimental39 studies concerning the dihydride facets in the SiNWs of different sizes. Bond lengths measured in the core layers (Sicore
Sicore) and at the surface (Si
H) of the equilibrated samples are equal to 2.441 ( 0.098 Å and 1.533 ( 0.007 Å, respectively. Thus, the predicted Sicore
Sicore bond lengths lie within 8% of the experimental values of 2.32
2.35 Å for bulk silicon.66,67 The reported values of Si
Si (2.393 ( 0.490 Å) and Si
H (1.523 ( 0.003 Å) bond distances obtained by calculations using semiempirical65,68,69 or first-principle27,28,63,70
72 techniques are close to our results. We should mention that the measured Si
Si distances in the surface layer are slightly shorter than those in the core layers. A similar shortening of the bonds was also found in the theoretical studies on H-passivated silicon nanowires.27,73 Such shorter Si
Si bond distances could be caused by the steric hindrance due to the hydrogen atoms. Surface layer is known to affect the stiffness and electronic properties of the SiNWs.29,65,74,75 Therefore, the contribution of the surface layer to the total volume is calculated and included in Table 1 as surface-to-volume ratio. If this ratio is large, the mechanical properties of the wire will be remarkably influenced by the surface layer. Consequently, one can expect a decrease of the surface effect from B1 to B3 (for the octahedral type of cross sections) and from A1 to A3 (for samples with tetrahedral cross sections). 3.2. Surface Layer under Tension. To show how the changes in the overall wire shape affect the behavior of the surface layer under tension, two pairs of wires of similar widths (such as A2, B1 and A3, B2) are discussed in this subsection. It is obvious that the facets with Si-monohydrides in the B-type wires are expected to have different response to tension than those in the A-type wires. To clearly understand the behavior of the surface layers under tension, geometrical parameters such as selected distances and dihedral and valence angles are measured within the surface layer at different facets. To ensure the specificity of such structural changes only for the surface layer (H and Si atoms), they are compared to the data collected from the core of the wire. 3.2.1. H-Passivation Layer. One of the ways to characterize the surface layer under tension is to show the changes of its H 3 3 3 H distances (see Figure 2a). We provide the snapshots taken for A3 and B2 wires in Figure 2a for an in-depth understanding. The top and side view snapshots depict the motion of H atoms relative to each other and to the {110} facets, respectively. For the sake of clarity, two typical Si
H bonds that are considered for measuring 12285

dx.doi.org/10.1021/jp201948g |J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C

ARTICLE

Figure 2. Structural characteristics of the surface layer of studied nanowires plotted for the applied strain. All the structural parameters such as distances, angles, and dihedral angles provided here are average values. (a) Distance between two H atoms on {100} and {110} facets monitored during tension. Related Si
H bonds are highlighted in red on each snapshot to show the direction of the Si
H bond used in calculations. (b) Position of the Si
H bond relative to the surface of the {110} facet characterized by the angle between the Si
H bond and the average plane formed by the nearest three Si atoms of the surface layer (solid lines). The H
Si
Si
H dihedral angle measured on the surface of the {100} facet characterizes the motion of two Si
H groups relative to each other on the surface of the {100} facet (dashed lines). (c) Distance between two Si atoms highlighted in red on snapshots in picture (a) that represent the motion of the Si part of the surface layer (solid lines). Changes of similar distances measured in core layers during tension are plotted by dashed lines. 12286

dx.doi.org/10.1021/jp201948g |J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C the distances between H atoms are higlighted in red here. Details of how we measured the H 3 3 3 H distances at different facets are provided in Figure S1 (see Supporting Information). From comparison of the starting points of the curves in Figure 2a, we notice that the values of H 3 3 3 H distances in the equilibrated samples of the same shape are always larger for the wires with a larger surface-to-volume ratio. This could be explained by the fact that the H atoms on the facets of smaller wires have less neighbors than larger wires, and thus, they have more space to adjust their positions. From the shape and location of the curves, one can recognize not only the facet-type dependent but also the shape-type dependent character of the motion of H atoms under tension. By comparison of the H 3 3 3 H distances of similar width A- and B-type wires of the {110} facets, we see a change in pattern from a continuous decrease in the A-type wires to an appearance of a slight increase at large strain values in the Btype wires. To get a 3D picture of the motion of H-passivation layer, the direction of motion of corresponding Si
H bonds is evaluated. For the H atoms at the {110} facets, we calculate the angle characterizing the position of the Si
H bond relative to the {110} facet itsef. Three nearest Si atoms in the surface layer are used to define the {110} facet. An increase of these angles (see solid lines in Figure 2b) reveals how the SiH groups respond to the transverse compression originated from the stretching of Si layers. High initial values of such angles in the well-relaxed {110} facets of A-type wires increase further by a maximum of 15 (A3 sample, solid black line) and 17 (A2 sample, solid dark cyan line). The values of similar angles calculated for the equilibrated B-type wires are lower than for the A-type wires due to the constraints caused by two neighboring {100} facets. Calculated angles in the B-type wires increase to the maximum of 20 (B1 sample, solid dark yellow line) and 24 (B2 sample, solid blue line) during stretching. To understand the movements of H atoms at the {100} facets of the B-type wires, we calculate the dihedral angle H
Si
Si
H (see Figure 2b, dashed lines) that is from Si
H of the two nearest SiH2 groups. The calculated dihedral angles are also plotted in Figure 2b. A continuous increase of these dihedral angles from slightly negative values (
10 and
5) to their maximum of ∼45 corresponds to the rotation of SiH2 groups. Such a rotation at the edges of the {100} facets leads to an elongation of H 3 3 3 H distances on the {110} facets of the octahedral samples after reaching the strain values of 0.10
0.20. 3.2.2. Si Atoms of the Surface Layer. Since the rotation of SiH2 groups has been detected during the stretching of wires, the behavior of Si atoms in the surface layer cannot be correctly guessed from the changes of related H 3 3 3 H distances. Therefore, we have calculated the Si 3 3 3 Si distances (solid lines in Figure 2c) for the SiH and SiH2 groups whose H 3 3 3 H distances have been discussed in the previous section. Similar Si 3 3 3 Si distances within the core of the wire (dashed lines in Figure 2c) have also been calculated to show the difference between the inner and outer layers. Regardless of the shape and size, the variation of all Sicore 3 3 3 Sicore distances for a particular facet can be described by the bunch of closely positioned dashed curves. Much more interesting information may be collected from the distribution of curves for the Si 3 3 3 Si distances measured at the surface. For example, it seems that a large contribution of Simonohydrides to the total surface results in more concave shape of the {110} facets. On the other hand, the predominance of Sidihydrides on the surface causes the swelling of the {110} facets.

ARTICLE

Figure 3. Snapshots of SiNW cross sections taken before ultimate strength was reached. We also provide the schematic description of motion of Si atoms within each cross section drawn based on the data from Figures S2 and S3 (Supporting Information).

Therefore, the {110} facet of A3 is the most concave, while the same facet in samples A2, B1, and B2 becomes more and more protuberant (see the structural representation in Figure 3). In contrast to the {110} facet, the Si 3 3 3 Si distances on the {100} facets increase with increase of strain. Since the direction in which these distances are measured is not strictly perpendicular to the [001] direction, we cannot use these distances as a direct evidence of the downward or upward movement of {100} facets in the B1 and B2 wires. However, the location of the red and magenta solid curves with respect to their core analogues (dashed lines) suggests the wire widthrelated (or surface dominated) amplification of Si 3 3 3 Si distance variation. For further understanding of the motion of Si atoms in the core and surface layers, some additional geometrical parameters have been measured and analyzed (see detailed explanations and Figures S2 and S3 in the Supporting Information). From these measurements, one can see that the presence of H atoms on the surface of the wires adds some strain in the direction perpendicular to the longitudinal one. To balance such strains, certain distortions occur depending on the overall shape of the wire and the type of facets. In the schematic part of Figure 3, red arrows represent the transverse motions of Si atoms in the core and surface layers during the stretching along the [001] direction. The size of the red arrows corresponds to the amplitude of motions of Si atoms in the core and surface layers. 12287

dx.doi.org/10.1021/jp201948g |J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C

ARTICLE

Table 2. Values of Young’s Modulus (E, GPa), Poisson’s Ratio (υ), Valence Angle Averaged Deformation before Rupture (VAAD, degrees), Ultimate Tensile Strength (GPa) and Strain, and Breaking Strain Calculated for Æ001æ Si Nanowires of Different Sizes wire code

E, GPa

υ

VAAD before rupture, degrees

ultimate tensile strength, GPa

ultimate tensile strain

A1

32.000 ( 1.145

0.014 ( 0.002







A2

50.517 ( 0.486

0.096 ( 0.002

18.78

15.9

0.365

0.480

A3

62.797 ( 0.380

0.116 ( 0.002

14.04

14.88

0.271

0.303

B1 B2

38.735 ( 0.638 64.531 ( 0.316

0.128 ( 0.002 0.195 ( 0.002

22.43 17.79

15.0 17.63

0.346 0.309

0.387 ∼0.35

B3

88.047 ( 0.266

0.205 ( 0.001







In all samples, the most intensive motion of Si atoms is registered across the {110} facets in both core and surface layers. In the A-type wires, the transverse compression happened mainly due to the surface layer of Si atoms. At the same time, the structure of {100} facets allows some relaxation during tension in the B-type wires. For example, two additional lines of Si atoms in the center of the {100} facet of the sample B2 are not affected by neighboring {110} facets, and therefore, move in the direction opposite to the direction of motion of Si atoms at the edges of the {100} facet. In contrast to the tetrahedral wires, both core and surface layers of the {110} facets are compressed in the octahedral wires. From the distribution and size of red arrows given in the schematic part of Figure 3, one can also see a relationship between the wire width and the degree of transverse contraction during tension. 3.3. Mechanism of Fracture. The crystalline structure of the silicon wire, which was cut along the [001] direction, does not have any Si
Si bond aligned with or being perpendicular to the direction of stretching. Therefore, no bond can be taken as the first and the most affected due to the tension. However, there are many valence angles that are located in the planes perpendicular to the longitudinal axis. These valence angles are significantly affected during the earlier period of stretching. When the valence angles reach their limit of increase or decrease, the bonds start to stretch and break; then the visible signs of fracture appeared. The valence angles that are located on the planes parallel to the [001] direction and perpendicular to the {110} facets have been measured in the core layers of A2, A3, B1, and B2 samples. These valence angles change about 15
22 before the fractures are seen. The data provided in Table 2 reveal the significant changes of valence angles for the small width wires having octahedral cross sections. The wire A1, which has the smallest width, shows an extra twisting of the structure. Here, the twisting becomes more powerful than the normal fluctuations caused by temperature. The combined effect of twisting and stretching results in extremely fast rupture of the wire A1 after the strain reached the value of 0.15. Two factors would lead to such an outcome. First, during the stretching, only the difference in distances between the pair of Si atoms at the edges of the wire is regulated, and no other constraint is applied. Second, the small width of the A1 wire results in a few points along the longitudinal direction, where the parts of the structure are connected by only one Si atom. Hence, the wire A1 has the freedom to rotate at those points during the stretching. Since the core structure in large width wires is similar to the bulk crystalline structure, there are no signs of significant twisting. Among all studied samples, the B3 wire should possess the most rigid structure of the core layers. We could simulate only the initial period of stretching because of the large size of the B3 system. Consequently, we cannot measure its valence angles before failure.




complete fracture strain










To show the details of structural changes up to the complete fracture, several snapshots for the A2, A3, B1, and B2 wires have been collected in Figure 4. Compared to the initial equilibrated structures (at strain = 0), the stretching causes a significant elongation along the [001] direction. In addition, it results in compression of all the wires (see the second snapshot for each wire) in the direction perpendicular to the longitudinal one. A cleavage mechanism is observed for all the wires except sample A2 wherein some Si atoms start to shear and form a new type of bonding after fracture is initiated. Starting from the third snapshot (see specific values of strain in Figure 4), the samples undergo successive bond breaking. We can only speculate about the preferred ways of rupture in nanowires due to the limited number of samples. For example, we observe that the B-type samples tend to break the bonds along the “>” line and form two almost complementary pieces after fracture. Moreover, the way in which the constraints are applied results in a slight shift of the relative positions of two pieces of the B2 sample just before the strain became equal to 0.347. After this event, the forces are applied along the line that has some deflection from the longitudinal axis. Therefore, we assume that the complete rupture of the sample B2 would have happened at the strain value of ∼0.35. This kind of situation was reported experimentally.19 We do not see any similarity in the fracture of A-type wires. Upon strecthing of the A3 sample, breaking of multiple bonds initiates propagation of two parallel cracks perpendicular to the [001] direction (strain = 0.293). This process results in local distortion of the structure and formation of a cavity in place of the first crack and a complete cleavage of all bonds along the second crack. Consequently, the broken pieces are not complementary. The A2 sample exhibits the most unusual mechanism of fracture. Initially, bonds start to break along the “>” line at the strain values of 0.376 and 0.378 and further formation of a strand-like connection between two parts of the wire (strain = 0.389 and 0.4). Finally, after the strand broke completely, the Si atom at one of its ends rebinds with neighboring atoms and forms triangular shape. At the same time, a loop of five SiH groups is formed at the surface of another piece of the wire. It is obvious that final pieces of the A2 wire are not complementary to each other. Interestingly, the fracture never occurred exactly at the center of the wire. We encounter four different mechanisms that suggest the presence of statistical uncertainty due to large availability of bonds. Silicon nanowires break at the strain values of 0.303
0.48 (see Table 2). The only available experiment related to the Æ001æ SiNW has been done for the wire that was partially passivated with the oxide layer, and its diameter is ∼6 nm. The experimental study reported the fracture at the strain value of 0.27 ( 0.01.19 Our study indicates that the complete fracture in larger width wires occurs earlier than in smaller width wires. In agreement with our observation, the fracture of larger width wire in experiment 12288

dx.doi.org/10.1021/jp201948g |J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C

ARTICLE

Figure 4. Snapshots of the samples A2, A3, B1, and B2 at different stages under tension. Si-atoms are marked in yellow and red, while H atoms are shown in gray. Red Si atoms and arrows represent place and direction of applied forces.

occurred earlier than the fracture of smaller width systems studied in this work. 3.4. Mechanical Properties. The relationship between stress and strain has been monitored during the simulation and presented in Figure 5. Orange and violet curves are incomplete because of the large size of the B3 system and the presence of excessive twisting in the smallest A1 sample. The stress
strain curves given in black (A3 wire) and dark yellow (B1 wire) colors start as linearly growing curves similar to

the elastic region of brittle materials. The subsequent plastic region is short, and then the curves fall abruptly due to the rupture process. Whenever the system has fewer connection points between its parts, there is a chance of tilt angle development between two pieces of the breaking wire because of the imperfection in the way of applying constraints. In such a situation, the uniaxial stretching could not be done perfectly. This problem has been noticed in the B2 wire, in which presumable rupture should occur at the 12289

dx.doi.org/10.1021/jp201948g |J. Phys. Chem. C 2011, 115, 12283–12292

The Journal of Physical Chemistry C

Figure 5. Stress
strain curves in the process of stretching. The  sign indicates complete fracture.

strain value of ∼0.35 and the stress value of ∼4 GPa. Furthermore, the blue curve starts to fluctuate around the stress value of 4 GPa for a while, and the B2 system completely broke at the strain value of 0.45. The stretching at strain values from 0.35 to 0.45 cannot be considered as the ideal stretching along the [001] direction. During this period, the right portion of the B2 sample (see Figure 4) departed and rolled over the left portion that is still attached to it by one Si
Si bond. A slightly different situation has happened with the A2 sample. The effect of shortcomings of the way of applying constraints is seen in the later period of the rupture of A2. We have also noticed that the dark cyan curve starts to fluctuate around the stress value of 4 GPa (see Figure 5). The fluctuating part of the stress
strain curve is related to the shear of atoms near the crack and the formation of new bonds to stabilize the breaking of A2 wire. In comparison with the blue curve (B2 wire), the dark cyan curve (A2 wire) has more fluctuations because of the reorganization of Si structure. Qualitatively, Figure 5 suggests that toughness and ductility increase when the sample width decreases. At the same time, the degree of stiffness and strength is high in wider samples. The quantitative measurements of mechanical properties such as Young’s modulus and Poisson’s ratio are calculated for SiNWs in the region of small strains (