Introduction to strain calculation

Here we introduce the theoretical background of the strain and stress calculation in nextnano++. At first we will describe the definition of a strain tensor \(\varepsilon\) and stress tensor \(\sigma\) and then describe the basis of strain tensor calculation in nextnano++. A strain tensor is used to calculate the shifts and splittings of band-edge energies and piezoelectric charges.

The detailed explanation for the syntax in strain{ } is here: strain{ }.

Table of contents

Strain tensor \(\varepsilon\)

The calculation of strain effects in nextnano++ is based on linear continuum elasticity theory, in which a crystal can be described by a field of material points with coordinates \(\mathbf{x}\). A distortion of the crystal shifts any point to a new position \(\mathbf{x}'=\mathbf{x}'(\mathbf{x})\). A field of displacement vectors \(\mathbf{u}\) is defined as the devision between the new position and the original position:

\[\mathbf{u}(\mathbf{x}):=\mathbf{x}'(\mathbf{x})-\mathbf{x}\]

Figure 2.2.1 The field of displacement vector \(\mathbf{u}\) at \(\mathbf{x}\). This is the vector along which the point that was at the position \(\mathbf{x}\) moved through the displacement.

A strain tensor \(\varepsilon\) is defined using this displacement vector:

\[\varepsilon_{ij}:=\frac{1}{2}\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right] ; \quad (i,j=1,2,3)\]

Strain is dimensionless. The diagonal elements of this strain tensor \(\varepsilon_{ii}\) represents the length changes per unit length in \(x_i\)-direction as described in Figure 2.2.2.

Figure 2.2.2 Deformation of a dilatable string in an unstrained (top) and strained state (bottom). We can see the diagonal element \(\varepsilon_{ii}=\frac{\partial u_i}{\partial x_i}\) represents the length changes per unit length in \(x_i\)-direction.

The off-diagonal elements \(\varepsilon_{ij(i\neq j)}\) arise due to shear deformations of the crystal. Figure 2.2.3 shows the deformation of an infinitesimal rectangle in \(x_1x_2\) plane. We can see \(\frac{\partial u_2}{\partial x_1}=\frac{u_2(x_1+\Delta x_1,x_2)-u_2(x_1,x_2)}{\Delta x}=\sin\alpha\simeq \alpha\) and \(\frac{\partial u_1}{\partial x_2}=\frac{u_1(x_1,x_2+\Delta x_2)-u_1(x_1,x_2)}{\Delta x_2}=\sin\beta\simeq \beta\). In these angle changes, \(\frac{\alpha-\beta}{2}\) corresponds to a pure solid-body rotation and \(\frac{\alpha+\beta}{2}=\frac{1}{2}\left[\frac{\partial u_2}{\partial x_1}+\frac{\partial u_1}{\partial x_2}\right]=\varepsilon_{12}\) measures the shear strain.

Figure 2.2.3 Deformation of an infinitesimal rectangle in a strained state.

By definition strain tensor \(\varepsilon\) is symmetric (i.e. \(\varepsilon_{ij}=\varepsilon_{ji}\)) so the number of components that must be specified is actually 6. Voigt notation is the useful convention in which these 6 independent components are written in form of a 6\(\times\)1 matrix for short. This notation reads:

\[11\rightarrow 1, 22\rightarrow 2, 33\rightarrow 3, 23 \rightarrow 4, 31\rightarrow5, 12\rightarrow6\]

and

\[\begin{split}\begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \varepsilon_{3} \\ \varepsilon_{4} \\ \varepsilon_{5} \\ \varepsilon_{6} \end{bmatrix} = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{23} \\ 2\varepsilon_{13} \\ 2\varepsilon_{12} \end{bmatrix}\end{split}\]

Stress tensor \(\sigma\)

A stress tensor component \(\sigma_{ij}\) represents the force towards \(x_j\)-direction acting on infinitsimal area that is perpendicular to \(x_{i}\)-direction. Its unit is the same with pressure (\(\mathrm{[Pa]}=\mathrm{[N/m^2]}\)).

Figure 2.2.4 The components of stress tensor \(\sigma\).

In linear approximation, this stress tensor is related to the strain tensor \(\varepsilon\) by means of Hook’s law:

\[\sigma_{ij} = \sum_{kl}C_{ijkl}\varepsilon_{kl}\]

where \(C_{ijkl}\) is the component of eleasticity stiffness tensor, which is the forth-order tensor comprising \(3^4=81\) components. It’s dimension is the same with stress tensor components and defined as \(\mathrm{[GPa]}\) in nextnano++. In Voigt notation, \(C\) is the form of a 6\(\times\)6 matrix by putting \(C_{ijkl}=C_{mn}\ (i,j,k,l=1,2,3,\ m,n=1,...6)\). Then the Hook’s law reads

\[\begin{split}\begin{bmatrix} \sigma_{1} \\ \sigma_{2} \\ \sigma_{3} \\ \sigma_{4} \\ \sigma_{5} \\ \sigma_{6} \end{bmatrix} = \begin{bmatrix} C_{11}&C_{12}&C_{13}&C_{14}&C_{15}&C_{16} \\ C_{21}&C_{22}&C_{23}&C_{24}&C_{25}&C_{26} \\ C_{31}&C_{32}&C_{33}&C_{34}&C_{35}&C_{36} \\ C_{41}&C_{42}&C_{43}&C_{44}&C_{45}&C_{46} \\ C_{51}&C_{52}&C_{53}&C_{54}&C_{55}&C_{56} \\ C_{61}&C_{62}&C_{63}&C_{64}&C_{65}&C_{66} \end{bmatrix} \begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \varepsilon_{3} \\ \varepsilon_{4} \\ \varepsilon_{5} \\ \varepsilon_{6} \end{bmatrix}\end{split}\]

For many crystal structures with high symmetry, many of these coefficients are 0 and some are related to others. The elasticity tensor of zincblende and wurtzite crystals are given by

\[\begin{split}C_{zb} = \begin{bmatrix} C_{11}&C_{12}&C_{12}& & & \\ C_{12}&C_{11}&C_{12}& & & \\ C_{12}&C_{12}&C_{11}& & & \\ & & &C_{44}& & \\ & & & &C_{44}& \\ & & & & &C_{44} \end{bmatrix}\end{split}\]

\[\begin{split}C_{wz} = \begin{bmatrix} C_{11}&C_{12}&C_{13}& & & \\ C_{12}&C_{11}&C_{13}& & & \\ C_{13}&C_{13}&C_{33}& & & \\ & & &C_{44}& & \\ & & & &C_{44}& \\ & & & & &C_{66} \end{bmatrix}\end{split}\]

with \(C_{66}=\frac{1}{2}[C_{11}-C_{22}]\) in wurtzite.

These constants are defined in database_nnp.in. You can also overwrite these values in your input file.

• For zinc-blend materials, for example:

  database{
      binary_zb{
          name = GaAs
          valence = III_V
  
          elastic_consts{
              c11 = 122.1                  # [GPa] elastic constants
              c12 = 56.6                   # 1 * 1011 dyn/cm2 = 10 GPa  ->  12.21 * 1011 dyn/cm2 = 122.1 GPa
              c44 = 60.0                   # The elastic constants are needed for the calculation of the strain in heterostructures.
          }
      }
  }
  

• For wurtzite materials, for example:

  database{
      binary_zb{
          name = GaN
          valence = III_V
  
          elastic_consts{
              c11 = 390                    # [GPa] elastic constants
              c12 = 145                    # 1 * 1011 dyn/cm2 = 10 GPa  ->  39.0 * 1011 dyn/cm2 = 390 GPa
              c13 = 106                    #
              c33 = 398                    #
              c44 = 105                    # The elastic constants are needed for the calculation of the strain in heterostructures.
          }
      }
  }
  

Strain and stress calculation

Next we will describe how the strain tensor \(\varepsilon\) and stress tensor \(\sigma\) are determined in general. Then the two types of calculation implemented in nextnano++ are introduced briefly.

In general

The principle of conservation of linear momentum results in the following equations of stress tensor components for \(i=1,2,3\):

\[\sum_{j=1}^3\frac{\partial{\sigma_{ji}}}{\partial x_j}+f_i=0\]

where \(\mathbf{f}\) is the body force such as gravity. When the boundary conditions are specified, the field of displacement vector \(\mathbf{u}\), by which the stress tensor components \(\sigma_{ij}\) are eventually written, is determined according to these simultaneous differential equations. Then the strain tensor \(\varepsilon\) and stress tensor \(\sigma\) are also determined from \(\mathbf{u}\).

Note

The principle of conservation of angular momentum, on the other hand, results in the symmetricity of stress tensor: \(\sigma_{ij}=\sigma_{ji}\)

The field of displacement vector which satisfies the above balance equations and boundary conditions also minimizes the total potential energy \(U+V_E\) where \(U\) is the elastic strain energy and \(V_E\) is the potential energy associated with the body force \(\mathbf{f}\). This is so called minimum total potential energy principle.

In the linear approximation regime, the elastic energy stored in the whole body is:

\[U=\frac{1}{2}\int_VC_{ijkl}\varepsilon_{ij}\varepsilon_{kl}\ dV\]

When the body force \(\mathbf{f}\) is assumed to be zero throughout the system, solving the above differential equations is equivalent to find the strain tensor that minimizes this elastic energy \(U\).

In nextnano++

There are two kinds of calculation of strain, pseudomorphic_strain{ } and minimized_strain{ }, in nextnano++. In both of implementations pseudomorphic layer is assumed as the boundary condition between the substrate and the layer grown on this substrate. The substrate is assumed to be so thick that the in-plane lattice constants of the layer is matched to that of substrate. Also, the body force \(\mathbf{f}\) is assumed to be 0 throughout the structure.

In this assumption, the analytic expressions for strain tensor that satisfies the aforementioned stress balance equations (i.e. that minimizes the elastic energy) can be found for 1D structures. This analytic solution is implemented on pseudomorphic_strain{ }. This feature also works in 2D or 3D but the user must be sure that the model makes sense from a physical point of view (i.e. the 2D/3D structure should consist of different layers along the growth direction whereas the layers must be homogenous along the two perpendicular directions).

On the other hand, minimized_strain{ } calculates the strain tensor by minimizing the elastic energy mentioned before. This can also be used for 1D simulations. In this case, the results will be equivalent to the analytical model pseudomorphic_strain{ }.

The detailed explanation for the syntax in strain{ } is here: strain{ }. Please refer to [AndlauerPhD2009] for more details about these topics.

Last update: nn/nn/nnnn