Atomistic modeling of the fracture toughness of silicon and silicon-silicon interfaces

Abstract

Fracture modeling at the atomic scale is currently an intense area of research because the crack propagation process depends strongly on the description of interatomic interactions. Here we present first the mode-I plane strain quasi-static fracture toughness of single crystal silicon, along four orientations, as obtained using molecular statics simulations with nine empirical potentials. The “best” potential is determined by comparing the fracture toughness and trapping range of the simulations with available experimental data and with results calculated using first-principles molecular dynamics. The choice is buttressed by its ability to predict the effective toughness and propagation direction of a crack subjected to mode-II loading. The best performing potential is then used to investigate the fracture toughness and the role of bond trapping for a crack along the boundary of two silicon crystals belonging to two different tilt families.

Appendices

Appendix A: Elasticity matrices

The elastic stiffness tensor of a (010)[001] oriented Si SC, that has cubic symmetry, can be expressed in matrix form using the Voigt notation as

$$\begin{aligned} \mathbf{C}^{(0)}= \left[ \begin{array}{cccccc} C_{11} &{} C_{12} &{} C_{12} &{} 0 &{} 0 &{} 0\\ C_{12} &{} C_{11} &{} C_{12} &{} 0 &{} 0 &{} 0\\ C_{12} &{} C_{12} &{} C_{11} &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} C_{44} &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} C_{44} &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} C_{44}\end{array} \right] , \end{aligned}$$

(A.1)

where $C_{11}$, $C_{12}$, and $C_{44}$ are the only three independent elastic constants, listed in Table 1 for different interatomic potentials. Next, the forth order stiffness tensor for an arbitrary crystal orientation $(\mathbf{n}_2)[\mathbf{n}_3]$, labeled in accordance with the crack system notation and defined by the rotation matrix ${{\varvec{\varOmega }}}$, can be obtained in the following way.

$$\begin{aligned} C_{ijkl}=\varOmega _{ip}\varOmega _{jq}\varOmega _{kr}\varOmega _{st}C^{(0)}_{ijkl}, \end{aligned}$$

(A.2)

where Einstein summation rule is used. The matrix ${{\varvec{\varOmega }}}$, that rotates the initial (010)[001] crystal system around axis $\mathbf{m}$ by angle $\theta _0$ so that it orients the crystal as $(\mathbf{n}_2)[\mathbf{n}_3]$ (i.e. it has $\mathbf{n}_2$ crystallographic orientation along y axis and $\mathbf{n}_3$ along z axis), writes as

$$\begin{aligned} \varOmega _{ij}=\cos \theta _0\delta _{ij}+\sin \theta _0 \epsilon _{ikj}m_k+(1-\cos \theta _0)m_im_j, \end{aligned}$$

(A.3)

where $\delta$ is the Kronecker delta and $\epsilon$ is the Levi-Civita symbol.

For a SC system we define the stiffness tensor $\mathbf{C}$ obtained by applying rotational transformation ${{\varvec{\varOmega }}}^{\mathbf{m}}_{\theta _0}$ shown in Eq. (A.3), where axis $\mathbf{m}$ and angle $\theta _0$ are dictated by the considered crystal orientation, to the initial $\mathbf{C}^{(0)}$ given in Eq. (A.1). For the tilt boundaries we also find the stiffnesses $\mathbf{C}^{(1)}$ and $\mathbf{C}^{(2)}$ obtained by additional rotation about z axis by angles $\theta _1$ and $\theta _2$ of a SC with stiffness $\mathbf{C}={{\varvec{\varOmega }}}^{\mathbf{m}}_{\theta _0}\mathbf{C}^{(0)}$; this transformation is denoted by ${{\varvec{\varOmega }}}^{z}_{\theta _{\beta }}$ where $\beta =1$ and 2 correspond to the solutions for $\theta \in (0,\pi )$ and $\theta \in (-\pi ,0)$ regions, respectively, see Fig. 1. Equivalently, $\mathbf{C}^{(\beta )}$ can be computed by applying the total rotation operation being a product of two as ${{\varvec{\varOmega }}}^{(\beta )}={{\varvec{\varOmega }}}^{z}_{\theta _{\beta }}{{\varvec{\varOmega }}}^{\mathbf{m}}_{\theta _0}$ to $\mathbf{C}^{(0)}$ according to Eq. (A.2). For example, a (110)[001] SC system has $\mathbf{m}=(0,0,1)^{{\mathsf {T}}}$, $\theta _0=\pi /4$; $(110)[1\bar{1}0]$ has $\mathbf{m}=(-0.357,-0.863,0.357)^{\mathsf{T}}$, $\theta _0=1.718$; $(111)[1\bar{1}0]$ has $\mathbf{m}=(-0.642,-0.762,0.085)^{\mathsf{T}}$, $\theta _0=1.578$; $(111)[11\bar{2}]$ has $\mathbf{m}=(-0.367,-0.887,-0.282)^{\mathsf{T}}$, $\theta _0=2.909$.

Appendix B: Linear elastic solution for a bimaterial interface crack

This section briefly reviews the linear elastic solution obtained by Qu and Bassani (1989) and Bassani and Qu (1989) and includes the key expression needed to implement the used boundary conditions. In a rectangular coordinate system $x_i$ ($i=1,2,3$; $x_1$, $x_2$, $x_3$ correspond x, y, z coordinate system used in the paper), the infinitesimal strain tensor, defined as $\varepsilon _{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})$, where comma refers to the partial derivative, $\mathbf{u}$ is the displacement, is related to the Cauchy stress tensor as $\sigma _{ij}=C_{ijkl}\varepsilon _{kl}$. Then the equilibrium equation writes as

$$\begin{aligned} C_{ijkl}u_{k,lj} =0. \end{aligned}$$

(B.1)

In case of the two-dimensional plain strain problem considered here $\varepsilon _{i3}=0$, $u_i=u_i(x_1,x_2)$. It allows to search for a solution in the form of $u_i=a_i f(z)$, where f is an arbitrary function of complex variable $z=x_1+px_2$. Substituting it into Eq. (B.1) leads to the following eigenvalue problem

$$\begin{aligned}{}[\mathbf{Q}+p(\mathbf{R}+\mathbf{R}^{\mathsf{T}})+p^2\mathbf{T}]\mathbf{a}=0, \end{aligned}$$

(B.2)

that yield six roots for p and corresponding eigenvectors a. Here are $\mathbf{Q}$, $\mathbf{R}$, $\mathbf{T}$ are the elasticity matrices defined as $Q_{ik}=C_{i1k1}$, $R_{ik}=C_{i1k2}$, $T_{ik}=C_{i2k2}$. Since $\mathbf{C}$ is positive definite, there are three pairs of complex conjugate roots, so that ${\text {Im}}(p_{\alpha })>0$, $p_{\alpha +3}=\bar{p}_{\alpha }$ (overbar denotes the complex conjugate) for $\alpha =1,2,3$. Next, additional matrices are introduced as $\mathbf{P}={\text {diag}}[p_1,p_2,p_3]=\langle p_{\alpha } \rangle$, $\mathbf{A}=[\mathbf{a}_1,\mathbf{a}_2,\mathbf{a}_3]$, $\mathbf{B}=\mathbf{R}^{\mathsf{T}}\mathbf{A}+\mathbf{T}\mathbf{A}\mathbf{P}$.

Moving to the interface crack problem, we first identify the matrices introduced above for a SC of a given orientation with stiffness $\mathbf{C}$ before tilt is applied, i.e. for perfect structure without the interface. Then, the properties of the system with a tilt boundary are computed as

$$\begin{aligned}&\mathbf{A}_{\beta }={{\varvec{\varOmega }}}^{x_3}_{\theta _{\beta }} \mathbf{A},\, \mathbf{B}_{\beta }={{\varvec{\varOmega }}}^{x_3}_{\theta _{\beta }} \mathbf{B},\, \mathbf{P}_{\beta }=\langle p^{(\beta )}_{\alpha } \rangle , \nonumber \\&\text {where} \qquad p_{\alpha }^{(\beta )}=\dfrac{p_{\alpha } \cos \theta _{\beta }+\sin \theta _{\beta }}{-p_{\alpha }\sin \theta _{\beta }+\cos \theta _{\beta }}, \end{aligned}$$

(B.3)

$\beta =1$ and 2 correspond to the solutions for materials having different tilt angles, i.e. for $y>0$ and $y<0$ regions, respectively. Solving for the asymptotic stress and displacement fields produced at the interface crack tip, Bassani and Qu (1989) introduced matrices $\mathbf{W}$ and $\mathbf{D}$ as

$$\begin{aligned} \mathbf{W}= & {} -{\text {Re}}(\mathbf{A}_1\mathbf{B}_1^{-1}-\bar{\mathbf{A}}_2\bar{\mathbf{B}}_2^{-1}),\nonumber \\ \mathbf{D}= & {} -{\text {Im}}(\mathbf{A}_1\mathbf{B}_1^{-1}-\bar{\mathbf{A}}_2\bar{\mathbf{B}}_2^{-1}), \end{aligned}$$

(B.4)

which reflect continuity of displacements and tractions across the bimaterial interface. Qu and Bassani (1989) also proved that if $\mathbf{W}=\mathbf{0}$ the crack tip fields are not oscillatory, i.e. solution for the stress field has the standard inverse square root singularity.

For the particular case of uniform tractions applied along the crack faces, that can be represented by the stress intensity vector $\mathbf{k} = (K_{\mathrm{II}},K_{\mathrm{I}},K_{\mathrm{III}})^{\mathsf{T}}$, for bimaterial satisfying $\mathbf{W}=\mathbf{0}$ the linear elastic solution for the stress being function of radial distance r and polar angle $\theta$ derived by Bassani and Qu (1989) takes the following form

$$\begin{aligned}&(\sigma _{12},\sigma _{22},\sigma _{32})^{\mathsf{T}}=\dfrac{{\varvec{\theta }}^{(\beta )}(\theta )\mathbf{k}}{\sqrt{2\pi r}}, \qquad \nonumber \\&(\sigma _{11},\sigma _{21},\sigma _{31})^{\mathsf{T}}=-\dfrac{\tilde{\varvec{\theta }}^{(\beta )}(\theta )\mathbf{k}}{\sqrt{2\pi r}}, \qquad \nonumber \\&{\varvec{\theta }}^{(1)}(\theta )\!=\!{\text {Re}}(\bar{\mathbf{B}}_1\bar{{{\varvec{\Lambda }}}}_1\bar{\mathbf{B}}_1^{-1}), \, \tilde{\varvec{\theta }}^{(1)}(\theta )\!=\!{\text {Re}}(\bar{\mathbf{B}}_1\bar{{{\varvec{\Lambda }}}}_1\bar{\mathbf{P}}_1\bar{\mathbf{B}}_1^{-1}),\nonumber \\&{\varvec{\theta }}^{(2)}(\theta )\!=\!{\text {Re}}({\mathbf{B}}_2{{\varvec{\Lambda }}}_2{\mathbf{B}}_2^{-1}), \, \tilde{\varvec{\theta }}^{(2)}(\theta )\!=\!{\text {Re}}({\mathbf{B}}_2{{\varvec{\Lambda }}}_2\mathbf{P}_2{\mathbf{B}}_2^{-1}),\nonumber \\&{{\varvec{\Lambda }}}_{\beta }=\langle {(\cos \theta +{p}_{\alpha }^{(\beta )}\sin \theta )^{-1/2}}\rangle . \end{aligned}$$

(B.5)

The asymptotic behavior of the displacements can be deduced from the full finite crack solution given in Eqs. (3.14a,b) in Bassani and Qu (1989) and it writes as

$$\begin{aligned}&{} \mathbf{u}(r,\theta )=\sqrt{\dfrac{2r}{\pi }}{\varvec{\theta }_\mathbf{u}}^{(\beta )}(\theta )\mathbf{k},\nonumber \\&{\varvec{\theta }_\mathbf{u}}^{(1)}(\theta )={\text {Re}}(\bar{\mathbf{A}}_1\langle {(\cos \theta +\bar{p}_{\alpha }^{(1)}\sin \theta )^{1/2}}\rangle \bar{\mathbf{B}}_1^{-1}),\nonumber \\&{\varvec{\theta }_\mathbf{u}}^{(2)}(\theta )={\text {Re}}({\mathbf{A}}_2\langle {(\cos \theta +p_{\alpha }^{(2)}\sin \theta )^{1/2}}\rangle {\mathbf{B}}_2^{-1}). \end{aligned}$$

(B.6)

The crack opening displacement can be calculated as

$$\begin{aligned} {{\varvec{\varDelta }}}(x_1)=\sqrt{2x_1/\pi }\mathbf{D}\mathbf{k} \, \text { for } \, x_1<0. \end{aligned}$$

(B.7)

Besides, according to Wu (1990), the energy release rate becomes

$$\begin{aligned} G_{c}=\dfrac{1}{4}{} \mathbf{k}^{\mathsf{T}}\mathbf{D}\mathbf{k}. \end{aligned}$$

(B.8)

Although this approach may seem complicated and unnecessary for the problems considered here, since the solution for a uniformly loaded interface crack of $\mathbf{W}=\mathbf{0}$ bimaterial is identical to the linear elastic solution for one anisotropic media obtained by Sih et al. (1965) with corresponding elastic properties for each region, it provides a more rigorous framework and allows to justify the solution type (oscillatory versus non-oscillatory) by computing the matrix $\mathbf{W}$.