Shear Instability in Nanoporous Si Joo-Hyoung Lee* School of Materials Science and Engineering, Gwangju Institute of Science and Technology, Gwangju, Republic of Korea ABSTRACT: Elastic properties of nanoporous Si (np-Si), which is composed of bulk Si containing ordered, nanometersized cylindrical pores, are investigated based on ﬁrst-principles density functional theory calculations. By separately varying the pore size and spacing, it is demonstrated that the elastic stiﬀness of np-Si under the shear strain perpendicular to the pore axis turns negative when the volume fraction of pores becomes greater than a critical value. The total energy calculations reveal that the negative values in the stiﬀness originate from the enhanced strain energy, which leads to signiﬁcant rotation in bonds near the pore surface. Moreover, the high sensitivity of the elastic stiﬀness to shear induces a structural transformation in np-Si from tetragonal (D2d) to orthorhombic (C2v) phase, which makes it necessary to properly take the eﬀect of external strain due to substrates or electrical leads into account in np-Si-based applications. KEYWORDS: Nanoporous Si, density functional theory, shear instability, phase transformation

W

compressed elastometric matrix with a square array of circular holes undergoes a pattern transformation to mutually orthogonal ellipses, which is caused by a reversible elastic instability in the porous structure.15,16 Bertoldi et al. further demonstrated that the strain-induced pattern transformation leads to an auxetic behavior (negative Poisson’s ratio),17 for which the Poisson’s ratio can be tuned by adjusting the aspect ratio of an alternating pattern of elliptical voids.18 On the theoretical side, Duan et al. and Ouyang et al. employed composite cylinder assemblage (CCA) model19 and generalized self-consistent method (GSCM) to show that the elastic constants of nanoporous materials are determined by the cylindrical pore size and porosity (ϕ, volume fraction of pores).20,21 Depending on the pore size and ϕ, it was also shown that the elastic constants can exceed those of the nonporous materials, which implies that nanoporous structures become stiﬀer than the parent materials. While all these ﬁndings contributed to elucidating the elastic responses of the materials containing cylindrical pores, however, either the pores are purely macroscopic (μm to mm in diameter) or results are based on continuum mechanics approach and attempts that take the atomistic nature of nanoporous structures into account have been lacking. In the present work, to gain an insight into the elastic properties of np-Si we calculate its elastic stiﬀness (C44 and C66) by employing ﬁrst-principles electronic structure calculations, which have been successfully applied to describe mechanical properties of such systems as metals,22−24 semiconductors,25,26 and nanostructures.27,28 Through separately varying the pore

ith rapid advances in fabricating sophisticated nanostructures, nanoporous materials have been receiving a great amount of attention in such diverse areas as ion exchange,1 gas storage,2,3 catalyst,4,5 molecular sieving,6,7 and water desalination applications8 due to the large volume fraction of pores and high surface areas. Among various nanoporous structures, those with an ordered arrangement of cylindrical pores have recently attracted much interest because of their high energy conversion eﬃciency in thermoelectric (TE) applications.9−11 Unlike the poor TE properties of the bulk phase, Si containing nanocylindrical pores (nanoporous Si, or np-Si) was demonstrated to display a sharply improved TE performance compared to the bulk at room temperature,9−12 which is due to the preserved electronic transport arising from the structural regularity of np-Si and increased phonon scattering. More interestingly, the TE conversion eﬃciency of np-Si can be further enhanced by applying external strain. Recent ﬁrstprinciples calculations have shown13 that the TE ﬁgure of merit at room temperature (ZT300 K) of strained np-Si is increased by as large as 70% from the unstrained value, which is caused by the modiﬁcation of the electronic structure due to the removal of degeneracy in the conduction band under strain. This ﬁnding is truly surprising in that ZT300 K of np-Si under strain becomes as high as that of Bi2Te3,14 one of the most widely used TE materials, thus making np-Si particularly attractive for solidstate cooling applications. However, it should be noted that all these studies are focused only on the electronic properties of np-Si, and another important characteristic of materials, mechanical properties, has yet to be clariﬁed for the nanoporous geometry. In regard to mechanical properties of materials containing cylindrical pores, Mullin et al. observed that a uniaxially © 2014 American Chemical Society

Received: May 13, 2014 Revised: August 13, 2014 Published: August 25, 2014 5081

dx.doi.org/10.1021/nl501772j | Nano Lett. 2014, 14, 5081−5084

Nano Letters

Letter

Figure 1. (a) Structure of np-Si with H-passivated pore surface, and (b) structural parameters for computing elastic constants. The pore axis is in zdirection.

diameter and spacing, it is demonstrated that both C44 and C66 approach the bulk value as the pore spacing becomes larger, whereas C66 signiﬁcantly softens for the increased pore diameter and eventually turns negative, which implies that the strained np-Si exhibits a structural instability and undergoes a permanent phase transformation. Total energy calculations further reveal that the instability arises from the enhanced strain energy stored in bonds, which leads to large bond rotations near the pore surface. The structure of np-Si employed in the present work is presented in Figure 1a, where the pore axis is oriented along [001] direction and the surface is H-passivated to remove dangling bonds. To examine the dependence of the elastic stiﬀness on the structural variables of pore diameter (dp) and spacing (ds) (Figure 1b), we consider 10 diﬀerent porous structures (S1−S10) as is listed in Table 1; from S1 to S5, ds is

In the case of np-Si with tetragonal symmetry, the relevant matrix elements required to calculate the stiﬀness constants are listed in Table 2.29 To obtain Cij’s, the total energy diﬀerences Table 2. Elements of the Strain Matrix D for Calculating Cija

ds (nm)

dp (nm)

ϕ

e0s (eV)

es (eV)

S1 S2 S3 S4 S5 S6 S7 S8 S9 S10

0.63 1.17 1.72 2.26 2.80 1.0 1.0 1.0 1.0 1.0

1.0 1.0 1.0 1.0 1.0 0.63 1.17 1.72 2.26 2.80

0.30 0.17 0.11 0.07 0.05 0.12 0.23 0.31 0.38 0.43

0.034 0.070 0.112 0.173 0.241 0.042 0.065 0.092 0.118 0.142

0.003 0.022 0.053 0.090 0.137 0.018 0.019 0.016 −0.149 −0.148

d1

d2

d3

d4

d5

d6

ΔE/V0

D1 D2

1 (1 + δ2)1/2

1 (1 + δ2)1/2

1 + δ2 1

δ 0

δ 0

0 δ

4C44δ2 2C66δ2

a ΔE is the diﬀerence in the total energy with and without strain, and V0 the volume of an unstrained unit cell, respectively.

between strained and unstrained structures are computed with four diﬀerent δ values (±0.005, ±0.01) and ﬁtted with a quadratic function. All calculations are carried out based on density functional theory (DFT) as implemented with Vienna Ab initio Simulation Package (VASP).31 The Kohn−Sham wave functions are expanded by employing plane waves with 350 eV cutoﬀ together with ultrasoft pseudopotential,32 and the exchange-correlation energy is treated within local density approximation parametrized by Perdew and Zunger.33 To perform the Brillouin zone integration, 24 irreducible k-points are considered with a Gaussian smearing parameter of 0.05 eV, respectively. In calculating the total energy of the strained structure, the atomic positions are relaxed while the unit cell is ﬁxed as deformed. We begin with the elastic properties of bulk Si, for which C44 (= C66) is computed to be 77 GPa, which is in excellent agreement with experiment.34 On the other hand, in the case of np-Si, the equivalence between C44 and C66 is lifted because the presence of pores reduces the crystal symmetry from the cubic Oh to tetragonal D2d phase. As is seen from Figure 2a, C44 starts at 36.8 GPa for S1 and gradually increases to as large as 66.6 GPa when (dp, ds) = (1.0, 2.80) nm because the material approaches the bulk phase. It is clear from Figure 2b that C66 with dp ﬁxed exhibits a similar behavior as C44, but its variation is more rapid than that of C44; while with (dp, ds) = (1.0, 0.63) nm C66 is as small as 5.3 GPa, which is less than 7% of the bulk value, it quickly grows to 57.5 GPa for ds = 2.80 nm. On the other hand, when ds is ﬁxed, C44 decreases from 53.5 to 35.5 GPa when the structure varies from S6 to S8, as is clear from Figure 2a, and the decrease is even more pronounced in case of C66; while C66 is 41.1 GPa for S6, it rapidly reduces to 14.0 GPa when dp = 1.71 nm. It is interesting to see that the dependence of C44 and C66 on the structural variables can be understood with the GSCM. As is seen from Figure 3, the computed elastic

Table 1. ds, dp, Porosity (ϕ), and the Strain Energy before (e0s ) and after (es) Relaxation under the Transverse Shear with δ = 0.005 for Each np-Si Structure structure no.

strain

varied with dp ﬁxed at 1.0 nm, and vice versa for the others. To compute C44 and C66, a small deformation is applied to the unit cell according to29,30

R′ = RD where R′ and R are 3 × 3 matrices composed of the lattice vectors with and without strain, respectively. D is the deformation matrix containing strain components, which is written as ⎛ d1 d6 d5 ⎞ ⎜ ⎟ D = ⎜ d6 d 2 d4 ⎟ ⎜⎜ ⎟⎟ ⎝ d5 d4 d3 ⎠ 5082

dx.doi.org/10.1021/nl501772j | Nano Lett. 2014, 14, 5081−5084

Nano Letters

Letter

strain remains positive for all δ values, indicating that the structure is stable against the strain. Around δ = 0, however, ΔE is found to become particularly ﬂat in that ΔE ≃ 15 and 33 meV for S9 and S10 for δ = ±0.005, respectively, whereas the corresponding values are as large as 267 and 345 meV in case of δ = ±0.01. The ﬂatness makes the higher order terms become important in expanding ΔE, which implies that deformed S9 and S10 lie outside the linear elastic regime. This is in sharp contrast with the structures from S1 to S8, where the energy diﬀerence is well described with a quadratic function of δ. Under the D2 strain (referred to as “transverse shear”), on the other hand, np-Si exhibits an instability for large pores, which is clear from the negative values of ΔE in Figure 4b. As is shown in Figure 4b, when |δ| ≤ δc = 0.0015 (0.001) for S9 (S10), |ΔE| is less than 5 meV, suggesting only a weak instability in this δ range. At δc, however, ΔE shows an abrupt drop (Δd) and then approaches a constant value (Δ0) as |δ| increases; for both S9 and S10 cases, Δd and Δ0 are found to be about 125 and 150 meV, respectively. The instability of np-Si with high porosity under the transverse shear arises from the rotation of the bond near the pore surface. In Figure 5, the changes in the bond angles of S9 as a function of the strain size are presented. It is found that among the angles between Si bonds, those at Si2−Si1−Si3, Si2−Si1−Si4, and Si3−Si1−Si5 bonds (Figure 5a) are the most susceptible to the shear. As is seen from Figure 5b, while the angle at Si2−Si1−Si3 and Si2−Si1−Si4 bonds are 107.3° and 110.4°, respectively, at equilibrium, they abruptly increase to 113.4° and 115.1° at δ = 0.0015, respectively, and continue to grow with δ. The angle at Si3−Si1−Si5 bond also shows a sharp change as δ varies but unlike the other two cases it exhibits a gradual decrease after a sudden reduction from the equilibrium value by 6° at δc. Similar changes are induced by the strain to the bond containing H atoms at the pore surface as is shown in Figure 5 (c). The bond angles at Si1−Si3−H1 and Si1−Si4−H3 increase by about 5° and 4° when δ = δc, respectively, whereas the angles become smaller by about 5° (Si1−Si3−H2) and 7° (Si1−Si4−H4) for the other H-containing bonds. It should be noted that similar variations in the bond angles are induced by the transverse shear for S10 for which more bonds are aﬀected due to the increased number of atoms. To understand the origin of the bond rotation, we calculate the strain energy before and after atomic relaxation. In Table 1, the initial (i.e., unrelaxed) strain energy (e0s ) with δ = 0.005 is listed for the nanoporous structures. As expected, e0s increases with the size of the unit cell; when dp is ﬁxed, e0s of S5 becomes as large as that of S1 by a factor of 7. On the other hand, in cases of ﬁxed ds the increase in e0s is slightly more than by a factor of 3 as dp varies from 0.63 to 2.80 nm. These initial strain energies are, however, substantially reduced upon the atomic relaxation as is seen from es in Table 1; in ﬁxed dp cases, the reduction is 91% for S1, gradually decreases as the cell size grows and becomes 43% for S5. In contrast, the reduction in the initial strain energy exhibits a linear increase with the unit cell when ds is ﬁxed; while e0s is reduced by 57% for S6 when the atomic positions are relaxed, the reduction becomes as high as 83% for S8. This implies that more initial strain energy is stored per bond with higher porosity, which is completely released through rotating bonds when ϕ reaches a critical value and thus induces a permanent structural deformation in S9 and S10. In summary, we have investigated the elastic properties of nanoporous Si by employing ﬁrst-principles electronic structure

Figure 2. (a) C44 and (b) C66 of np-Si. Circles (squares) are for the cases of ﬁxed dp (ds).

Figure 3. Calculated elastic stiﬀness of np-Si (dots) together with ﬁtting function from GSCM: (a) C44 and (b) C66.

stiﬀnesses are well ﬁtted with the results from the continuum approach; for C44, the ﬁtting function takes a form of35 0 C44 = C44

6(1 − ϕ) + 6σ(1 − ϕ) + κ(4 − 3ϕ) + κσ(2 + 3ϕ) 2[3 + 2ϕ + σ(3 − 2ϕ) + κ(2 + ϕ) + κσ(1 − ϕ)]

where κ = α1/(C044dp) and σ = α1/(C044dp) and for C66 C66 = β2

0 (1 − ϕ)d pC66 + (1 + ϕ)β1 0 (1 + ϕ)d pC66 + (1 − ϕ)β1

is employed,20 respectively. Here, C044 (= C066) is the elastic stiﬀness of bulk Si, and the ﬁtting parameters are obtained as (α1, α2) = (−8.68, 75.76) and (β1, β2) = (−43.88, 88.0) for C44 and C66, respectively. Interestingly, C66 appears to be less than zero when ϕ is large, as can be seen from Figure 3b. Indeed, C66 becomes negative if dp > 2 nm when ds is ﬁxed (S9 and S10). To understand the negative values in the elastic stiﬀness, we calculate the total energy diﬀerences (ΔE) between strained and unstrained cases of S9 and S10 and present this in Figure 4. As is seen from Figure 4a, ΔE of S9 and S10 under the D1

Figure 4. Total energy diﬀerence (ΔE) between strained and unstrained np-Si under (a) D1 and (b) D2 strains, respectively. Circles and squares are for S9 and S10, respectively. 5083

dx.doi.org/10.1021/nl501772j | Nano Lett. 2014, 14, 5081−5084

Nano Letters

Letter

Figure 5. (a) Atomic arrangement of S9 with atom indices near the pore surface; (b) angles between Si bonds, Si2−Si1−Si3 (circle), Si2−Si1−Si4 (square), and Si3−Si1−Si5 (up triangle); and (c) angles between the bonds containing H atoms, Si1−Si3−H1 (circles), Si1−Si3−H2 (square), Si1− Si4−H3 (up triangle), and Si1−Si4−H4 (down triangle). (11) Tang, J.; Wang, H.-T.; Lee, D. H.; Fardy, M.; Huo, Z.; Russel, T. P.; Yang, P. Nano Lett. 2010, 10, 4279. (12) Lee, J.-H.; Grossman, J. C.; Reed, J.; Galli, G. Appl. Phys. Lett. 2007, 91, 223110. (13) Lee, J.-H. Phys. Chem. Chem. Phys. 2014, 16, 2425. (14) Snyder, G. J.; Toberer, E. C. Nature Nater. 2008, 7, 105. (15) Mullin, T.; Deschanel, S.; Bertoldi, K.; Boyce, M. C. Phys. Rev. Lett. 2007, 99, 084301. (16) Willshaw, S.; Mullin, T. Soft Matter 2012, 8, 1747. (17) Bertoldi, K.; Reis, P. M.; Willshaw, S.; Mullin, T. Adv. Mater. 2010, 22, 361. (18) Taylor, M.; Francesconi, L.; Gerendás, M.; Shanian, A.; Carson, C.; Bertoldi, K. Adv. Mater. 2014, 26, 2365. (19) Hashin, Z.; Rosen, B. W. J. Appl. Mech. 1964, 31, 223. (20) Duan, H. L.; Wang, J.; Karihaloo, B. L.; Huang, Z. P. Acta Mater. 2006, 54, 2983. (21) Ouyang, G.; Yang, G.; Sun, C.; Zhu, W. Small 2008, 4, 1359. (22) Roundy, R.; Krenn, C. R.; Cohen, M. L.; J. W. Morris, J. Phys. Rev. Lett. 1999, 82, 2713. (23) Clatterbuck, D. M.; Krenn, C. R.; Cohen, M. L.; J. R. Morris, J. Phys. Rev. Lett. 2003, 91, 135501. (24) Kamran, S.; Chen, K.; Chen, L. Phys. Rev. B 2009, 79, 024106. (25) Roundy, R.; Cohen, M. L. Phys. Rev. B 2001, 64, 212103. (26) Zhang, R. F.; Veprek, S.; Argon, A. S. Phys. Rev. B 2008, 77, 172103. (27) Ertekin, E.; Chrzan, D. C. Phys. Rev. B 2005, 72, 045425. (28) Teich, D.; Fthenakis, Z. G.; Seifert, G.; Tománek, D. Phys. Rev. Lett. 2012, 109, 255501. (29) Mehl, M. J.; Osburn, J. E.; Papaconstantopoulos, D. A.; Klein, B. M. Phys. Rev. B 1990, 41, 10331. (30) Ravindran, P.; Fast, L.; Korzhavyi, P. A.; Johansson, B.; Wills, J.; Eriksson, O. J. Appl. Phys. 1998, 84, 4891. (31) Kresse, G.; Furthmueller, J. Phys. Rev. B 1996, 54, 11169. (32) Vanderbilt, D. Phys. Rev. B 1990, 41, 7892. (33) Perdew, J. P.; Zunger, A. Phys. Rev. B 1981, 23, 5048. (34) Hopcroft, M. A.; Nix, W. D.; Kenny, T. W. J. Microelectromech. Syst. 2010, 19, 229. (35) Zhang, W. X.; Wang, T. J. Appl. Phys. Lett. 2007, 90, 063104.

theory. Our calculations have demonstrated that the nanoporous structure becomes unstable against the shear strain perpendicular to the pore axis, which is evidenced from the negative C66 for the porosity greater than a critical value. The structural instability of np-Si is found to originate from the bond rotation due to the enhanced strain energy, which in turn leads to an irreversible transformation to a based-centered orthorhombic structure. It should be noted that the shear instability is also observed for other pore orientations, implying that it is a generic response of highly porous np-Si under shear strain. This observation makes it important to examine the possibility of structural deformation due to external strain exerted by substrates or electrical leads which are in contact with the porous structure. Because the physical properties of np-Si are likely to be modiﬁed due to such permanent distortion, our ﬁnding provides a valuable guide in designing relevant structures with optimal porosity for future np-Si-based applications.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing ﬁnancial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MSIP) No. 2013R1A1A2013412.

■

REFERENCES

(1) Plabst, M.; McCusker, L. B.; Bein, T. J. Am. Chem. Soc. 2009, 131, 18112. (2) Millward, A. R.; Yaghi, O. M. J. Am. Chem. Soc. 2005, 127, 17998. (3) Morris, R. E.; Wheatley, P. S. Angew. Chem., Int. Ed. 2008, 47, 4966. (4) Hu, Y.-S.; Guo, Y.-G.; Sigle, W.; Hore, S.; Balaya, P.; Maier, J. Nat. Mater. 2006, 5, 713. (5) Sharma, K. K.; Anan, A.; Buckley, R. P.; Ouellette, W.; Asefa, T. J. Am. Chem. Soc. 2008, 130, 218. (6) O’hern, S. C.; Boutilier, M. S.; Idrobo, J.-C.; Song, Y.; Kong, J.; Laoui, T.; Atieh, M.; Karnik, R. Nano Lett. 2014, 14, 1234. (7) Joshi, R. K.; Carbone, P.; Wang, F. C.; Kravets, V. G.; Su, Y.; Grigorieva, I. V.; Wu, H. A.; Geim, A. K.; Nair, R. R. Science 2014, 343, 752. (8) Cohen-Tanugi, D.; Grossman, J. C. Nano Lett. 2012, 12, 3602. (9) Lee, J.-H.; Galli, G. A.; Grossman, J. C. Nano Lett. 2008, 8, 3750. (10) Yu, J.-K.; Mitrovic, S.; Tham, D.; Varghese, J.; Heath, J. R. Nat. Nanotechnol. 2010, 5, 718. 5084

dx.doi.org/10.1021/nl501772j | Nano Lett. 2014, 14, 5081−5084